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the  computat ional  plane. 

Ratio  of  circumference  to  diameter  of  a  circle. 

Fluid  density. 

Similarity  parameter  in  Blasius'  series  of  Equation  16; 
and  a  grid  transformation  coefficient  defined  *>v  Equation 
A .  7  . 

Inner  and  outer  control  volume  surfaces. 

Grid  transformation  coefficient  defined  by  Equation  A.h. 
Element  of  the  viscous  stress  tensor. 

Element  of  the  turbulent  stress  tensor. 

Element  of  the  total  stress  tensor,  '  .  +  ... 

i  .1  t  ‘  J 

Transonic  small  disturbance  potential  doublet  fun,  tion. 
Transonic  small  disturbance  potential  vortex  iun.tiin. 

Acceleration  parameter  for  pressure  iteration  give-  ’  . 

Equation  4o  defined  as  2  .!“/((  i  4 

pb 

.  n  v  i  -  s  i  d  st  rear,  function. 


xv  i 


Svmbo 1 

Magnitude  of  the  local  flow  field  vorticity,  curl  vq  and 
a  local  acceleration  parameter. 

Acceleration  parameter  for  the  field  pressure  finite 
^  difference  equations. 

,  Acceleration  parameter  for  the  surface  pressure  iteration, 

pb 

w  Acceleration  parameter  for  the  u  and  v  velocity  component 

uv 

finite  difference  equations. 


Subscript 

b 

e 


f 

i 


j 


k 


Denotes  a  location  on  the  body  surface. 

Denotes  a  location  at  the  edge  of  the  boundary  layer  or 
wake  . 

Denotes  a  location  on  the  computational  far  field  boundary. 

Denotes  the  f  coordinate  location  in  the  computational 
plane  for  finite  difference  equations;  and  an  indicial 
notation  index  in  differential  equations. 

Denotes  the  r,  coordinate  location  in  the  computational 
plane  for  finite  difference  equations;  and  an  indicial 
notation  index  in  differential  equations. 

A  component  index  in  indicial  notation. 


max ,min 


Denote  the  maximum  and  minimum  values  of  a  variable. 


S  Denotes  a  value  at  the  laminar  separation  point  which 

defines  the  beginning  of  the  bubble. 


w 

Denote; 

■i  a  value  in  the 

wake 

region . 

x,xx 

Denote 

dif  f  erent iat ion 

with 

respee t 

to  X 

y  >  y  v 

Denote 

d i f  f  erent iat ion 

with 

respec t 

to  y 

T:,rv 

Denote 

d i f  f  erent iat ion 

with 

respect 

t  o  f 

r  r  f 

’  * 

Denote 

d i f  f  erent iat ion 

with 

respect 

to 

Denotes  a  vector  quantity. 


Supe r sc  r ip t 

*  Denotes  a  quantity  in  the  transformed  plane. 

Denotes  a  turbulence  time  averaged  quant i tv. 


x  v  i  i 


cript 

Denotes  a  turbulent  fluctuating  quantity. 

Denotes  a  temporary  nondimensional  variable  definition 

Iteration  number  in  the  successive-over-relaxation 
iteration  procedure. 

Pertaining  to  P  contours. 

Pertaining  to  £  contours. 

Derivative  operator. 

Substantial  derivative  operator  (D-/Dt  =  c* — / t  +  v*  . -) 
Increment . 

Summation . 

Del-operator  defined  as  i  +  im¬ 
partial  derivative  operator. 
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ABSTRACT 

Numerical  solutions  are  obtained  for  two-dimensional  incompressible 
turbulent  viscous  flow  over  airfoils  of  arbitrary  geometry.  An  algebraic 
eddy  viscosity  turbulence  model  based  on  Prandtl's  mixing  length  theory 
is  modified  for  separated  adverse  pressure  gradient  flows.  Finite  dif¬ 
ference  methods  for  solving  the  inviscid  stream  function  equation  and  the 
incompressible  laminar  N'avier-Stokes  equations  are  used.  A  finite  difference 
method  for  solving  the  Reynolds  averaged  incompressible  turbulent  two- 
dimensional  N’avier-Stokes  equations  is  employed. 

The  inviscid  stream  function  equation  and  the  Navier-Stokes  equations 
are  transformed  using  a  curvilinear  transformation.  A  body-fitted  coor¬ 
dinate  system  with  a  constant  coordinate  line  defining  the  airfoil  section 
surface  is  transformed  to  a  rectangular  coordinate  system  in  the  trans¬ 
formed  or  computational  plane.  An  elliptic  partial  differential  Poisson 
equation  for  each  coordinate  is  used  to  generate  the  coordinate  system  in 
the  physical  plane  for  arbitrary  airfoils. 

The  two-dimensional  time  dependent  Reynolds  averaged  incompressible 
Navier-Stokes  equations  in  the  primitive  variables  of  velocity  and 
pressure  and  a  Poisson  pressure  equation  are  numerically  solved.  Turbu¬ 
lence  is  modelled  with  an  adverse  pressure  gradient  eddy  viscosity  tech¬ 
nique.  An  implicit  finite  difference  method  is  used  to  solve  the  set  of 
transformed  partial  differential  equations.  The  system  of  linearized 
simultaneous  difference  equations,  at  each  time  step,  is  solved  using 
successive-over-relaxation  iteration.  Far  field  boundary  conditions  are 


I 


examined.  Solutions  for  a  NACA  0012  airfoil  at  angles  of  attack  varying 
from  five  to  11.5  degrees  at  a  chord  Reynolds  number  of  170,000  are 
obtained.  Velocity  profiles  near  the  airfoil  surface  and  surface  pres¬ 
sure  distributions  are  presented  and  compared  with  experimental  data. 

Lift  and  drag  coefficients  agree  well  with  experimental  values.  The 
computed  lift  coefficients  near  stall  are  within  5%  of  the  experimental 
measurements,  and  the  numerical  drag  coefficients  agree  within  ten  drag 
counts  in  the  region  of  maximum  lift  to  drag  ratio.  The  short  laminar 
separation  bubble  near  the  suction  pressure  peak  is  numerically  determined. 
The  variation  of  bubble  length  and  turbulent  transition  length  with  angle  of 
attack  are  similar  to  experimental  trends. 
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SECTION’  I 


INTRODUCTION 


Much  effort  has  been  expended  by  the  aeronautical  commur.i  t  v  in  d,  - 
termining  the  aerodynamic  characteristics  of  airfoils.  1. inear  method; 
are  extensively  used  in  design  work  for  studying  configurations  at  sr.a 1 1 
angles  of  attack  with  negligible  flow  separation.  Experimental  wind 
tunnel  investigations  are  used  to  determine  the  characteristics  near  stall 
where  separation  phenomena  become  important.  Recent  developments  in 
numerical  techniques  have  stimulated  research  on  another  approach,  namelv 
the  numerical  solution  of  the  Navier-Stokes  equations.  These  equations 
model  the  viscous  effects  which  contribute  to  airfoil  stall.  For  this 
reason,  numerical  Navier-Stokes  methods  offer  the  possibility  of  determin¬ 
ing  the  aerodynamic  characteristics  for  airfoils  experiencing  stall. 
Numerical  methods  can  also  complement  experimental  methods  by  efficiently 
extending  the  range  of  parameters  under  investigation.  Furthermore, 
numerical  methods  eliminate  model  support  interference  and  wall  interference 
effects  found  in  wind  tunnel  testing. 

The  purpose  of  this  investigation  is  to  develop  a  numerical  Navier- 
Stokes  method  that  will  accurately  determine  the  aerodynamic  character ist  ics 
of  incompressible  turbulent  viscous  flow  over  two-dimensional  airfoils  near 
stall. 

The  development  of  a  numerical  method  for  turbulent  flow  requires 
a  suitable  technique  for  distributing  points  throughout  the  flow  field 
and  a  model  which  describes  the  behavior  of  turbulence  within  spicific 
regions  of  the  flow  field.  A  survey  of  numerical  grid  generating  l eehn i  ;ue;. 
and  available  turbulence  models  is  presented.  The  quantitv  of  1  iter.i!i;ii 
concerned  with  the  Navier-Stokes  equat ions  is  extensive.  The  ret  r<  ,  .< 
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summary  of  the  literature,  which  describes  formulations  of  the1  Navier- 
Stokes  equations  and  their  numerical  solving  techniques  applied  to  flows 
over  airfoils,  is  given.  The  research  objectives  for  this  work  are 
then  discussed,  and  a  summary  of  the  remaining  sections  is  presented. 

Body-fitted  curvilinear  coordinate  systems  greatly  enhance  the 
application  of  numerical  methods  to  practical  boundary  value  problems 
involving  partial  differential  equations.  The  representation  of  a 
boundary  surface  as  a  coordinate  line  reduces  the  difficulties  associated 
with  numerically  specifying  boundary  conditions  by  interpolation  in  finite 
differences  methods.  In  the  physical  (x,yl  plane,  values  for  one  compu¬ 
tational  coordinate  t  are  specified  at  selected  points  on  both  the  body 
surface  and  the  outer  boundary.  Constant  values  for  the  other  coir.putaf.iona 
coordinate  r  are  specified  on  both  the  body  surface  and  the  outer  flow- 
boundary.  The  transformed  computational  (;,  ’)  plane  then  becomes  a 
rectangular  region  with  an  orthogonal  grid.  Winslow  (3)  and  Chu  {2) 
introduced  the  concept  for  two-dimensional  regions  interior  to  a  elost.l 
boundary.  Their  transformed  coordinates  are  solutions  of  Laplace's 
equation  in  the  physical  plane  and  define  a  triangular  mesh  system  in 
the  physical  plane.  Amsden  and  Hirt  (31  took  the  physical  coordinates 
to  he  solutions  of  a  modified  Laplace's  equation  in  the  transformed 
plane . 

Thompson,  Thames,  and  Mast  in  (-1,3,  M  generalized  the  method  for  the 
automatic  generation  of  hodv-fitted  coordinates  for  any  two-d imens i ona 1 . 
mul t i-connected  region.  They  also  introduced  the  use  of  forcing  iunction.- 
in  a  Poisson  equation  which  provides  nesh  control  in  regions  with  large 
gradients.  Hodge  (71  developed  an  automated  grid  line  attraction  method 
based  cn  houndary  layer  theory  which  determines  the  coefficients  in  his 


forcing.  function.  Hodge  assumed  a  Blasius  boundary  layer  profile  and 
distributed  his  grid  points  at  approximately  equal  velocity  increments 
in  the  boundary  layer.  Steger  and  Sorenson  (S)  introduced  auxiliary 
conditions  for  the  forcing  functions  which  provide  ancle  and  distance 
control  at  the  inner  boundary  surface.  The  angle  with  which  a  line 
intersects  the  body  surface  is  specified  by  a  function  ('),  and  the 
rate  of  change  of  arc  length  with  ’  on  a  :  line  at  tin-  bodv  surface  i  - 
prescribed  by  s.  (  0  ■  Sorenson  (9)  later  imposed  similar  conditions  on 
an  outer  computational  boundary.  These  auxiliary  conditions  are  used  to 
solve  for  coefficients  in  the  chosen  exponential  forcing  fund  ions .  'he 
geometric  conditions  hold  exactly  only  in  the  limit  as  ’  approaches  arc. 
Sorenson  reported  that  numerical  instabilities  occur  for  large  change.- 
in  the  coefficients  during  successive  iterations  and  for  boundaries  with 
sharp  corners.  He  implemented  a  limit  function  which  damped  the  change  of 
the  value  for  each  coefficient  over  one  iteration;  and  at  sharp  corners  he 
computed  average  values  of  each  coefficient  from  data  at  the  neighboring 
boundary  points.  Mast  in  and  Thompson  (101  have  also  extended  tin.  elliptic 
body-fitted  coordinate  generation  technique  to  three  dimensions  for  simple 
geomet  r ies  . 

Some  other  specific  grid  generation  techniques  which  use  properties  of 
elliptic  differential  equations  have  also  beer,  reported.  Meydcr  (111 
constructed  an  orthogonal  curvilinear  coordinate  system:  by  using  electric 
field  theory.  He  solved  the  potential  equation  twice  with  different  sets 
of  mixed  Diriehlet  and  Neumann  boundary  conditions  for  the  electric  potential 
and  electric  force  lines,  respectively.  These  solutions,  however,  were 
obtained  in  the  physical  plane  using  in  t  erpola  t  i  or,  and  became  tlu  eu  r  •  i  1  i:  ,  a  r 
coordinate  lines  in  the  physical  plane.  Tlu  coordinate  metrics  were  thin 
used  to  formulate  a  finite  difference  equation  which  was  als,  solved  in  the 
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phvsical  plane.  Conformal  mapping  techniques  which  employ  the 


'1  i.eojorbtn- 

Garick  (11)  trar.s:  Tiaation  have  been  examined  by  Ives  (13).  He  introduced 
the  use  of  fast  Fourier  transform  methods  and  developed  a  new  class  of 
transformat  ions  which  naps  the  flow  field  of  a  two-element  airfoil  onto 
the  region  between  two  concentric  circles.  Conformal  techniques  are  nut 
capable  of  being  extended  to  three-dimensional  geometries  (4).  Neither 
of  these  approaches  offers  a  convenient  means  for  mesh  control  in  regions 
of  large  flow  gradients. 

Recently,  a  geometric  grid  generation  technique  has  been  introduced  by 
(libeling,  Shamroth,  and  Eisenan  (14).  The  technique  parameterizes  (t)  the 
body  and  outer  boundary  surfaces  and  uses  a  stretching  function  R(r)  for 
mesh  attraction  along  constructed  lines  perpendicular  to  the  body  surface. 
I'nit  increments  for  ordered  pairs  (t,r)  generate  the  corresponding  computa¬ 
tional  plane.  Further  refinements  which  provide  angle  and  arclength 
variation  control  at  a  body  surface  have  subsequently  been  developed  by 
E iseman  (15,16,17)  . 

The  search  for  an  accurate  and  universal  turbulence  model  has  paralleled 
the  development  of  body-fitted  grid  generation  techniques.  In  principle, 
the  Navi er-Stokes  equations  completely  describe  the  turbulent  fluctuating 
fluid  motion.  The  required  mesh  resolution,  necessary  to  resolve  the  tur¬ 
bulent  eddies  with  varying  length  scales,  when  translated  into  com: uter 
resources  presently  make  this  approach  unfeasible.  Many  quantities  of  en¬ 
gineering  interest  in  turbulent  flows  involve  a  mean  value  taken  over  a  time 
interval.  The  time  interval  is  sufficiently  long  to  include  man"  fluctuations 
while  snail  compared  with  the  characteristic  time  of  the  mean  flow.  The 
Navi  er-Stokes  equations  can  then  be  re-formulated  using  these  r.t  an  flow 
variables.  This  Reynolds  averaging  procedure  introduces  the  near,  of  terms 
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involving  fluctuating  quantities. 


The  Reynolds  stress  components 


-i  u!u!  are  the  most  common  terms  of  this  tvpe.  In  order  to  solve  the 
i  J 

averaged  form  of  the  Navier-Stokes  equations,  "turbulence  closure"  must  be 
achieved  by  suitably  modelling  these  additional  terms.  This  approach  to 
turbulence  has  led  to  models  which  introduce  auxiliary  relationships  ranging 
from  algebraic  equations  to  several  partial  differential  equations.  These 
models  are  commonly  categorized  and  are  now  summarized  according  to  the 
number  of  additional  partial  differential  equations  which  comprise  the  model. 

The  algebraic  or  zero  equation  models  have  their  origins  in  Bouss inesr. ' s 
(18)  eddy  viscosity  hypothesis  and  Prandtl's  (19)  mixing  length  model.  The 
local  turbulent  stresses  are  assumed  proportional  to  the  local  mean  flow 
strain  rates  with  the  proportionality  constant  defined  as  an  equivalent 
or  eddy  viscosity.  The  eddy  viscosity  models  of  Cebeci-Smith  (20),  Xeller- 
Herring  (21),  and  Patankar-Spalding  (22)  represent  this  approach.  The 
boundary  layer  in  each  method  is  divided  into  an  inner  near  wall  layer  and 
an  outer  wake  layer  with  separate  expressions  for  the  eddy  viscosity  co¬ 
efficient.  In  the  Cebeci-Smith  and  Patanker-Spal ding  models  the  inner  mixing 
length  varies  as  a  linear  function  of  the  normal  distance  from  the  wall  modi¬ 
fied  with  a  Van  Driest  (23)  laminar  sublayer  correction.  The  Cebeci  and  Mellor 
models  both  base  the  outer  length  scale  on  the  displacement  thickness  while 
the  Patanker-Spalding  model  uses  the  boundary  layer  thickness  directly. 

These  models  have  been  successfully  extended  and  applied  (24)  to  a  wide  variety 
of  boundary  layer  flow  geometries  involving  compressibility,  heat  and  mass 
transfer,  and  curvature  effects.  In  addition,  the  mean  flow  models  are 


computationally  efficient.  Launder  and  Spalding  (25)  point  out,  however, 


that  these  models  predict  a  vanishing  eddy  viscosity  where  the  velocity 
gradient  is  zero  and  have  not  been  successful  for  large  separated  recirculat¬ 
ing  flows. 


5 


The  one  equation  turbulence  model  introduced  by  Prandtl  (26)  is  an 


extension  of  the  algebraic  technique.  In  this  approach  the  solution  of 

the  partial  differential  turbulent  kinetic  energy  equation  usually  provides 

a  local  velocity  scale  given  by  q"  =  u|u|  where  summation  is  implied. 

The  length  scale  L  is  prescribed  by  an  algebraic  expression  as  before.  Terms 

involving  fluctuating  quantities  other  than  q  must  still  be  modelled. 

Glushko  (27),  Mellor  and  Herring  (28),  and  Wolfshtein  (29)  model  the  gradient 

diffusion  tern  -((u1^  pi-  )  +  u'^_  (u*^  u'j  )/ 2)  with  a  gradient  of  q  given 

as  NT  •'  ’  (c'/2)  where  N_  is  a  specified  constant  or  function  and  ■  ,,  is 

Q  M  k  Q 

the  eddy  viscosity.  The  Reynolds  stress  is  related  to  the  mean  flow  as 
in  the  algebraic  approach  except  that  the  eddy  viscosity  is  assumed  propor¬ 
tional  to  qL.  Bradshaw,  et  al  (30)  model  the  gradient  diffusion  tern,  with 
an  expression  G  q'  where  G  is  an  empirical  constant  or  function  and 
is  a  velocity  characteristic  of  large  eddy  motions.  They  also  write  the 

Revnolds  stress  directly  as  a  function  of  q  in  the  xorm  7  ,.=A..  q'.  In 

tij  ij  1 

each  model  the  isotropic  dissipation  is  modelled  by  a  form  q  7L.  Nee  and 
Kovasznay  (31)  use  a  rate  equation  for  the  total  viscosity  N  =  +  e., 

in  place  of  the  kinetic  energy  equation.  They  formulated  expressions  for 
the  generation  and  dissipation  terms  involving  various  constants  and 
the  length  scale.  Launder  and  Spalding  (25)  observe  that  the  one  differential 
equation  models  require  a  moderate  increase  in  computer  resources  but  do 
not  in  general  provide  more  accurate  results  than  the  results  obtained  from 
algebraic  methods. 

Kolmogorov  (32)  introduced  the  somewhat  more  complex  type  of  turbulence 
model  which  uses  two  partial  differential  equations.  A  form  of  the  turbulent 
kinetic  energy  equation  provides  a  local  velocity  scale.  The  local  length 
scale  is  obtained  from  a  second  equation.  Ng  and  Spalding  (33)  formalized 


6 


this  approach  by  introducing  the  energy-length  equation  derived  by 
Rotta  (34).  They  also  used  the  Glushko  closure  model  in  the  turbulent 
energy  equation.  This  '.lass  of  two  equation  models  is  named  the  k  -  kL 
turbulence  model  where  the  eddy  viscosity  is  proportional  to  k‘  L.  Saffman 
(35)  has  used  a  transport  equation  for  the  mean  vorticity  •  together  with 
the  turbulent  kinetic  energy.  The  eddy  viscosity  becomes  proportional  to 
k/~  in  this  k  -  -  class  of  closure  models  for  turbulence.  The  local  length 

J 

scale  is  assumed  to  be  proportional  to  k "/ ■ •  Wilcox  and  Rubesin  (36)  haw 
modified  this  approach  for  compressible  flows  and  generalized  the  constit¬ 
utive  equation.  Jones  and  Launder  (37)  use  a  transport  equation  for  the  rat 
of  dissipation  of  turbulent  energy  e  along  with  the  turbulent  kinetic 
energy  equation.  Glushko  type  closure  models  are  assumed  for  the  terms  in¬ 
volving  fluctuating  quantities  in  both  equations.  This  k  -  class  of  two 

3/2 

equation  models  has  a  local  length  scale  proportional  to  k  "  with  the 
eddy  viscosity  proportional  to  k  /-;  .  The  two  equation  models  have  been 
applied  to  various  boundary  layer  and  free  shear  layer  flows  with  a  variety 
of  constants  and  closure  models.  The  two  equation  models  require  signifi¬ 
cant  increases  in  computational  resources  and  have  not  led  to  a  model  of 
universal  applicability  (24,25). 

The  search  for  a  more  general  turbulence  model  has  led  to  the  use  of 
the  transport  equations  for  the  Reynolds  stresses.  In  this  approach  the  edd 
viscosity  concept  is  discarded.  Closure  models,  however,  are  still  required 
for  the  terms  containing  fluctuating  quantities  other  than  the  Reynolds 
stresses.  Donaldson  (38)  introdu.  closure  models  which  express  the 
fluctuating  terms  as  functions  of  the  Reynolds  stress  and  chosen  length 
scales.  Hanjalic  and  Launder  (39)  introduce  closure  models  which  use 
the  Reynolds  stress  and  retain  the  equations  for  turbulent  kinetic  energy 
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Only  cases  with  the  simplified  two- 


k  and  turbulent  dissipation  w  . 
dernensional  boundary  layer  approximations  have  been  investigated.  The 
extensive  computational  resources  and  initial  state  of  development  of  these 
models  preclude  this  type  of  approach  from  practical  consideration  as  a 
turbulence  model  for  a  complex  flow  field  computation. 

Theoretical  turl  lienee  models  of  further  complexity  have  appeared. 

Kolovandin  and  Yotutin  (40)  introduced  a  statistical  theory  where  additional 
equations  are  obtained  for  other  correlations  of  the  fluctuating  quantities. 
Ferziger  (41)  took  a  meteorological  viewpoint  of  turbulent  flows  by  numeri¬ 
cally  simulating  large  scale  eddies  while  modeling  the  small  scale  structures 
with  an  eddy  viscosity  technique.  This  approach  is  a  first  step  toward 
numerically  solving  a  turbulent  flow  field  with  the  instantaneous  equations 
of  motion. 

The  review  of  available  turbulence  models  provides  the  basis  for  select¬ 
ing  a  suitable  approach  for  use  in  the  numerical  solution.  The  algebraic 
and  the  two  equation  eddy  viscosity  models  appear  to  be  realistic  choices. 

As  previously  reported  (25),  the  one  equation  techniques  yield  unimproved 
results  compared  with  algebraic  methods  and  at  additional  cost  in  computer 
resources.  The  two  equation  methods  require  the  solution  of  additional  partial 
differential  equations.  In  these  methods,  terms  involving  fluctuating  quanti¬ 
ties,  except  the  Reynolds  stresses,  must  still  be  modelled  using  additional 
coefficients.  For  these  reasons,  the  simpler  algebraic  eddy  viscosity  technique, 
which  requires  considerably  less  computer  resources,  is  employed.  If  the 
physical  phenomena  associated  with  separated  adverse  pressure  gradient  flows 
can  be  included  in  the  turbulence  model,  then  the  accurate  calculation  of  the 
aerodynamic  characteristics  may  be  accomplished  with  an  algebraic  technique. 

A  survey  of  the  previous  research  involving  numerical  solutions  of 
the  N'avier-Stokes  equations  for  flow  over  airfoils  establishes  the  pre- 
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diction  level  of  computational  methods  and  also  reviews  the  numerical 
algorithms.  The  grid  generation  technique  and  turbulence  model  employed, 
when  applicable,  are  also  included. 

Early  numerical  solutions  of  the  Navier-Stok.es  equations  for  flow  over 

airfoils  used  the  vorticity-stream  function  formulation  with  an  automated 

grid  generation  technique  (4).  Walker  (42)  applied  the  method  to  the  laminar 

flow  over  a  flat  plate  and  compared  the  numerical  solution  with  Blasius' 

(43)  solution.  Thames  (44)  used  body-fitted  coordinates  with  the  vorticitv- 

stream  function  approach  and  solved  the  Navier-Stokes  equations  for  various 

bodies  in  doubly  connected  regions.  He  obtained  solutions  for  the  flow  over 

/ 

airfoils  at  chord  Reynolds  numbers  less  than  10  .  Problems  mainly  attributed 
to  wall  vorticity  developed  for  solutions  of  airfoils  at  angle  of  attack. 

Mehta  and  Lavan  (45)  also  used  a  vo  "ticity-stream  function  method 
in  studying  the  laminar  starting  vortex  and  separation  bubbles  for  impul¬ 
sively  started  incompressible  laminar  flow  over  a  Joukowski  airfoil  at 

4 

Reynolds  numbers  less  than  10  .  They  used  three  point  backward  time 
differences  and  centered  spatial  differences.  The  computational  grid 
was  obtained  through  a  conformal  transformation  followed  by  a  radial 
stretching  transformation. 

Reddy  and  Thompson  (46)  applied  an  integro-differential ,  vorticity- 
velocity  field  method  for  the  solution  of  incompressible  flow  in  doubly 
connected  regions.  Backward  time,  centered  spatial  (BTCS)  differences 
were  applied  to  the  Navier-Stokes  equations.  The  difference  equations 
were  solved  using  successive-over-relaxation  (SOR)  iteration.  They 
also  employed  the  coordinate  mesh  attraction  technique  with  a  time  dependent 
expanding  mesh  system.  Symmetric  airfoils  at  zero  angle  of  attack 
with  a  Reynolds  number  less  than  10^  were  considered.  At  the  higher 
Reynolds  numbers,  the  calculation  of  surface  vorticitv  required  a  large  mm- 
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her  of  grid  points  on  the  surface  which  greatly  increased  the  required 
computer  time.  A  steady  state  solution  was  not  obtained. 

The  vorticitv-velocity  field  formulation  has  been  uppliid  by  Sankar 
and  U’u  (47)  to  the  case  of  incompressible  laminar  f 1 ov  about  an  oscillating 
airfoil.  They  used  a  12'  thick  Joukowski  airfoil  at  a  Reynolds  nut hi r  of 
'  1000.  Triangular  finite  elements  were  constructed  near  the  airfoil  surface 
and  a  rectangular  mesh  was  used  away  iron  the  airfoil.  A  conformal  trans¬ 
formation  was  used  to  map  the  airfoil  in  the  physical  plant-  onto  a  unit  circle 
in  the  computational  plane.  Sankar  and  Tassa  (48)  investigated  this  same 
problem  with  a  compressible  flow  formulation  of  the  Navier-Stokos  equation.-.. 
They  used  a  conformal  transformation  followed  by  an  algebraic  radial  1  ret  ch  ir. . 
transfer::. at  ion.  The  alternating-direction-implicit  (AFT)  finite  difference 
method  of  Briley  and  Mci'onald  (491  was  used  to  obtain  solutions  for  Reynolds 
numbers  less  than  10  and  a  Mach  number  of  0.2.  In  a  separate  research 
effort,  Sugavanam  and  Wu  (50)  attempted  to  use  a  two  equation  k-  turbu¬ 
lence  model  with  a  vorticitv-velocity  formulation.  A  conformal  transforma¬ 
tion  for  a  12  percent  thick  Joukowski  airfoil  was  used.  They  exper iexced 
difficulties  in  obtaining  a  converged  solution  even  after  eight  hours  of 
CYBEk  74  CPU  time  were  expended.  Variation.-,  for  both  lift  and  drag  of  the 
order  of  50  percent  occurred. 

The  use  of  the  primitive  variables  of  velocity  and  pressure  for 
incompressible  flow  was  introduced  by  Harlow  and  Void,  (51)  in  the 
explicit  forward  time,  centered  spatial  (FTC'S)  M.irk-and-Col  1  ()”./"»  method. 

They  included  a  Poisson  equation  for  pressure  whirl,  is  obtained  b;.  taking 
the  divergence  of  the  momentum  equation.  They  also  f cmr.d  that  the  velocity 
divergence  terms  in  the  momentum  equation--  wi  r.  r<  qui  red  for  the  pressure 
field  calculation.  Hirt  and  Harlow  (r>2)  further  dove  loped  tin  r.et  b-  -d . 

Hodge  (53)  considered  the  case  of  laminar  ittv  or.pressibK  viscou.  :  1 .  -w  1  or 


an  airfoil  at  angle  of  attack.  He  used  a  for"  ■  !  t'a  T! 
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body-fitted  grid  transformation  and  applied  an  implicit  BTCS  differencing 
method  to  the  Navier-Stokes  equations.  The  system  of  difference  equations 
was  solved  with  SOR  iteration.  Various  methods  of  calculating  the  pressure 
field  were  investigated.  Hodge  concluded  that  the  momentum  equations  should 
retain  the  velocity  divergence  terms  and  that  a  Poisson  pressure  equation 
should  be  used  to  satisfy  continuity.  Ghia,  Hankey,  and  Hodge  (34)  applied 
a  Poisson  pressure  equation  and  the  primitive  variables  form  of  the  Navier- 
Stokes  equations  to  study  incompressible  driven  flow  in  a  square  cavity 
for  Reynolds  numbers  under  1000.  They  used  an  alternating-direction-implicit 
(ADI)  finite  difference  technique  for  the  momentum  equations  and  SOR  iteration 
for  the  pressure  equation.  A  Neumann  boundary  condition  derived  from  Lie- 
normal  component  of  the  momentum  equations  was  employed  to  compute  the  wall 
pressure . 

The  Poisson  pressure  equation  has  been  further  examined  by  Chien 
(55)  for  internal  flows.  He  observed  that  the  calculation  of  pressure 
by  direct  integration  of  the  momentum  equations  can  be  inaccurate  when 
large  velocity  variations  are  present.  He  found  that  forms  of  the 
pressure  equation  were  more  suitable  for  computing  the  pressure  field 
near  boundaries.  Recently,  Hodge,  et  al  (7)  applied  an  implicit  backward 
time  finite  difference  method  with  both  upwind  and  centered  spatial 
differences  to  the  Navier-Stokes  equations  in  primitive  variables  form. 

A  Poisson  pressure  equation  was  used  to  obtain  a  solution  for  laminar 
flow  over  airfoils  at  angle  of  attack  with  a  Reynolds  number  of  order  10  . 

The  solution  contains  a  large  oscillating  separated  region  which  gives 
approximate  trends  in  lift  but  very  poor  agreement  with  available  drag  data 
(50)  . 

Parly  numerical  Navi  r-Stokes  investigations  of  turbulent  flow  over 
airfoils  considered  two-dimensional  transonic  viscous  flow  at  small  angle*. 

1  1 


of  attack.  Deiwert  (57)  applied  MacCormack  s  explicit  method  (58)  with 


a  50x38  rectangular  based  exponentially  stretched  mesh.  He  used  20  uni¬ 
formly  spaced  points  to  define  the  upper  airfoil  surface  and  imposed  the 
Neumann  boundary  layer  condition  dp/Sn  =  0  on  the  surface.  Steger  (59) 
used  the  implicit  Beam-Warming  method  (60)  with  a  71x33  grid  obtained  from 
a  modified  Thompson  transformation  (4).  He  calculated  the  solution  on  the 
branch  cut  with  a  linear  extrapolation  procedure  and  evaluated  surface  pres¬ 
sures  using  the  momentum  equations  for  the  direction  normal  to  the  surface. 
Walitt,  et  al  (61)  applied  the  first  order  method  of  Trulio  (62)  coupled  with 
a  boundary  layer  technique  used  near  the  airfoil  leading  edge.  They  used  a 
130x68  mesh  system  but  did  not  capture  the  full  .suction  pressure  on  the  upper 
surface  near  the  airfoil  nose.  Each  investigation  used  a  Cebeci-Smith  (/ 5) 
type  algebraic  eddy  viscosity  model  and  the  Reynolds  averaged  compress :M e 
Navier-Stokes  equations. 

Recently,  Shamroth  and  Gibeling  (64)  used  the  Br i ley-McDona 1 d 
(4 5)  implicit  finite  difference  formulation  with  a  constructive  type 
81x30  grid  system  to  compute  turbulent  flow  over  an  airfoil  at  ar.  an.-li 
of  attack  of  6  degrees  and  Reynolds  number  of  10^.  They  first  tried 
a  two  equation  k-c  turbulence  model  but  experienced  convergence  problems 
near  the  leading  edge  and  far  field  flow  regions.  They  reported  obtaining 
large  or  negative  values  for  the  turbulent  viscosity  in  essent ial 1 v 
laminar  regions.  The  turbulent  kinetic  energy  equation  with  an  algebra i cal ly 
prescribed  length  scale  was  then  used.  Major  discrepancies  occurred  for  the 
mean  pressure  distribution  solution  in  the  suction  peak  and  trailing  edge 
regions  of  the  airfoil  surface.  They  also  noted  that  the  construct ive 
grid  technique  caused  a  crossing  of  grid  lines  of  the  same  family  for 
highly  cambered  airfoils. 

The  examination  of  existing  numerical  solutions  of  the  Nnvier-btokes 
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equations  for  two-dimensional  flow  over  airfoils  reveals  Reynolds  number 
and  angle  of  attack  restrictions.  As  a  result,  numerical  Navier-Stokes 
methods  for  the  accurate  computation  of  the  flow  field  and  resulting 
aerodynamic  characteristics  for  an  airfoil  near  stall  are  unavailable.  The 
development  of  such  a  numerical  method  for  solving  incompressible  two-dimen 
sional  turbulent  flow  over  airfoil  sections  near  stall  does,  therefore, 
constitute  a  significant  contribution  in  computational  aerodynamics.  The 
fulfillment  of  this  research  goal  requires  the  formulation  of  an  adequate 
turbulence  model  for  use  in  both  leading  and  trailing  edge  separated  region 
The  far  field  boundary  conditions  for  incompressible  flows  must  also  he 
examined . 

The  selection  of  a  suitable  numerical  approach  for  incompressible- 
turbulent  flows  is  based  on  several  considerations.  An  implicit  technique 
is  preferred  because  of  the  required  small  grid  spacing  near  the  body  sur¬ 
face  necessary  to  resolve  viscous  stresses.  Stability  criteria  of  explicit 
methods  would  impose  a  small  time  step  restriction  as  a  result  o',  the  arid 
sice.  The  numerical  method  should  he  capable  of  predicting  laminar  flow 
for  Reynolds  numbers  approaching  turbulent  flow  conditions  since  regions  of 
laminar  flow  nay  occur.  The  use  of  primitive  variables  is  readily  extended 
to  three-  dimensions  and  provides  for  the  direct  calculation  of  the  prcs-ur. 
field.  These  criteria  art-  satisfied  by  an  implicit  finite  difieri-r.ee  pro¬ 
cedure  developed  b\  Hedge  (Vi,?)  which  use!,  primitive-  v.iria:  !e  and  since 
over  -  re  1  axa  t  i  on  iteration.  The  implicit  iirite-  difference-  met  h  i-  «  -  ;  lo¬ 
in  this  investigation. 

The  rer.ei  in  inv  sect  ions  of  t  h  i  s  research  e  f  t  o  r  t  d  i  s,  . .  s  the  te  hi  by.. 

•r.p  loved  in  achieving  the  defined  research  voa  ’.  .  !  •  Si.  tie:.  !  1  a  !  r-  .  • 

the  Thompson  numerical!-,  generated  hodv-t  itt  «••!  «  o.  rdinati  t  •  '  i  >h  ■  l  : 1 

The  t  h  id  of  attracting  gr  id  line-  the  1  -d  -urf.c,  a-  d  d.-trm  :t  ;• 
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surface  points  are  also  discussed.  The  governing  equations  for  incon- 
pressible  two-dimensional  turbulent  flow  in  primitive  variables  are 
presented  in  Section  III.  Boundary  and  initial  conditions  are  discussed. 

An  algebraic  eddy  viscosity  turbulence  model  formulated  for  adverse 
pressure  gradient  separated  flows  is  described.  An  algebraic  type  of 
turbulence  model  has  been  chosen  because  of  the  previous  computational 
successes  for  other  turbulent  flow  problems  and  the  definite  computational 
efficiency  obtained  by  this  approach.  The  method  for  evaluating  force 
coefficients  exerted  on  the  body  by  the  flow  field  is  then  discussed. 
numerical  techniques  used  in  obtaining  a  solution  are  presented  in  Sect ie:. 

IV.  First,  the  finite  difference  method  for  numerically  generatin'  tin 
grid  is  given.  The  N'avier-Stokes  equations  are  then  written  in  finite 
difference  form.  The  numerical  method  used  to  implement  the  boundary  and 
initial  conditions  is  described.  The  computational  procedure  for  introducing 
turbulence  is  next  discussed.  The  finite  difference  and  numerical  inte¬ 
gration  techniques  for  calculating  the  force  coefficients  are  then  presented . 
The  formal  discussion  of  the  numerical  solutions  and  comparisons  wit:  exp.  r i- 
mental  data  are  given  in  Section  V.  The  conclusions  derived  fro:-,  thi-  rest  ar.  1 
work  and  recommendations  tor  future  work  art  set  forth  in  Section  Vi. 


Partial  derivatives  are  transformed  using  the  chain  rule  and  invert 


transformation  relations.  For  a  sufficiently  dif fei entiahle  function  f 
of  x  and  y,  the  derivatives  transform  as  follows: 

f  =  (v  f  .  -  v  .  f  )/J  ^ ' 

X  '  r  i  '  ; 

fv  =  (“xTl  f  t  +  x  r  frW.I  (V) 

Higher  order  derivatives  are  similarly  derived. 

The  system  of  elliptic  coordinate  generating  equations  is  chosen  to 
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and  on  curve  C  . 
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where  y)  and  v)  are  prescribed  functions  and  V  and  '  are 

prescribed  constants.  The  other  two  sides  C  and  C,,  in  the  rectangular 
region  are  transformed  from  the  branch  cut  -  Cn.  The  functions  >:(  7,  ’1 
and  y(  r)  and  all  derivatives  are  continuous  across  this  cut.  These 
boundary  conditions  insure  that  the  body  and  the  far  field  boundaries  art 
defined  by  constant  ’  lines.  The  generalized  functions  P  and  Q  art  utiliri 
in  attracting  the  coordinate  lines  in  various  regions. 

The  rectangular  grid  system  in  the  transformed  plane  is  used  for  the 


comptit  at  ions .  The  mesh  generating  system,  liquations  f>a  and  (<!> ,  is  trans¬ 
formed  to  the  computnt ional  plane  and  becomes 
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near  the  body  surface  require s  a 

£  in«.'  mesh  spacing  in  this  region.  The  forcing  functions  F  and  y  are  use 
to  attract  ’  and  ’  lines,  respectively.  In  this  work  the  ’  line  spacing 
is  accom.p  !  i.-he  d  by  sp»  c  i  fying  the  boundary  condition  j  unction  ’,(y,  yl  . 


The  •  line  spacing  is  adjusted  through,  the  use  of  the  0  function  propose 
by  Th.or.p~on  (-'.l  and  modified  by  Hodge  (71 
K  (  ') 

OC  ,  •  1  =  -  '  A.  (  ~)  exp. -P(  I’--,  1  ( 1 

k=l  k 

where  are  the  amplification  factors,  D  is  the  damping  factor,  K  is  th. 
number  of  terms,  and  •  ,  =  ( k  -  1 V  •  .  The  amplification  and  damp ing  1  act 

requiring  that  the  tangential  ve loci  tv  u  he  ;i  linear 
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are  determined  hv 


The  second  part  of  the  boundary  layer  grid  line  attraction  technique 


involves  calculating  values  of  and  D  which  when  used  in  0  give  the 

desired  n  line  spacing,  hear  the  body  surface,  the  T  transf ormat icn 

equation,  Equation  6b,  can  be  written 

(’1  )  =  Q(  h  ’•  )  (20) 

yyb 

where  variations  with  respect  to  x  are  neglected.  If  the  expression  for 
the  forcing  function  Q  in  Equation  15  is  substituted  into  Equation  20, 
then  Equation  20  can  be  integrated  with  respect  to  r  to  obtain 
..  2  K 

(^7)  =  k=l  exp[-D  H  (r  ,  ^k)]  +  H(sgn(’‘k  -  ^ )  ,  0.) 

[1  -  exp  (-P  r  -  hk  )]';  (21) 

where  the  difference  function  H  and  sgn  function  are  defined  as  follows 
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if  r,  >  r 
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sgn(rk  -  r)  = 


1  if  r -  r,  >  0 
k  — 


-1  if 


<  0 


Now  evaluate  the  expression  for  (dh/dy)'  ,  Equation  21,  at  r  =  0  and 
0  =  r;  -  which  gives  the  approximate  expressions  written  in  Equation 
22  and  23,  respectively. 
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(23) 


The  derivative  in  Equation  22  can  be  evaluated  using  the  y  values  obtained 


previously.  The  derivative  in  Equation  23  can  be  specified  by  choosing  an 
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Ti  based  on  the  placement  of  the  outer  boundary.  Then  Equations  22 
and  23  can  be  combined  to  give  the  following  equation  for  the  damping 
factor  D: 


D  =  In  2(~-) 
dv 


(— ) 
=  0  a\' 


•  .  7  •  -  V-  -  r 

-  ’  max  }. 


(24) 


With  the  damping  factor  D  known,  the  analytic  expression  (Equation  21) 
can  be  evaluated  on  each  interval  of  successive  v  values  obtained  from 
the  Blasius  solution  along  with  the  last  Li  interval  near  the  outer  boun¬ 
dary.  This  procedure  gives  a  system  of  K  equations  for  the  K  unknowns 
A,  which  are  then  solved.  Therefore,  values  for  the  amplification  factors 
and  damping  factor  are  obtained  at  each  r  line  which  are  used  to  compute 
a  body  fitted  grid  system. 

C .  Body  Surface  Point  Distribution 

The  body  surface  can  be  defined  by  either  a  given  set  of  points  or 
by  an  analytical  expression.  The  body  surface  grid  points  can  then  be 
determined  by  either  interpolation  or  evaluation  of  the  analytical  expres¬ 
sion. 

The  NACA  0012  airfoil  section  is  used  in  this  investigation.  This 
airfoil  has  been  used  in  several  investigations  (65)  and  is  an  AGARD  (66) 
designated  test  airfoil  section.  The  NACA  four  digit  series  airfoils 
have  an  analytic  description  of  the  body  surface  given  by  (67) 

y  (x)  =  5t  ■  0. 2969x  ‘  -  0.  126x  -  0.3516.x2  +  0.2843x3  -  0. 1 01 5x“'/  (25) 

where  x  is  the  distance  along  the  chord  (nondimensionalized  by  chord),  y 
is  the  upper  surface  ordinate,  and  t  is  the  thickness  given  as  a  fraction 
of  the  chord.  For  the  NACA  0012  airfoil,  t  =  0.12.  The  thickness  function 
given  by  F.quation  25  does  not  close  the  body  at  the  trailing  edge  (x=l). 
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A  circular  arc  is  used  to  close  the  body.  The  arc  is  constructed  to  have 
a  center  on  the  chord  line,  to  have  points  of  tangency  with  the  upper  and 
lower  airfoil  surfaces  near  x  =  0.99,  and  to  intersect  the  chord  line  at 
x  =  1.0. 

The  distribution  of  grid  points  on  the  airfoil  surface  is  based  on 
both  resolving  streamwise  flow  field  gradients  and  defining  the  surface 
curvature.  These  criteria  require  the  clustering  of  points  in  the  nose 
and  trailing  edge  regions  of  the  airfoil.  The  following  analytic  trans¬ 
formation  is  used  to  obtain  the  desired  clustered  distribution  of  x 
abscissa  values 

A3 

tanh(- — r^1')  -  tanh(-  ~) 

x  =  •  - *= - -  (2b) 

tanh( — — )  -  tanh( — — ) 

where  z  represents  the  input  set  of  equally  incremented  values  (0  ^  z  <  1), 
x  represents  the  clustered  set  of  abscissae  (0  <  x  <  1),  and  Al,  A2,  and  A3 
are  arbitrary  constants.  The  constant  Al  determines  the  center  of  the 
hyperbolic  function  and  is  used  to  cluster  points  toward  one  end.  The  con¬ 
stant  A2  varies  the  slope  of  the  function  at  the  center  and  thus  determines 
the  distribution  near  mid  chord.  Tiie  exponent  A3  can  also  affect  the 
clustering  of  points  at  the  ends.  If  A3  >  1,  more  points  occur  near  x  =  0; 
while  if  A3  <  1,  points  are  clustered  more  at  x  =  1. 
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SECTION  III 


GOVERNING  EQUATIONS 


The  time  dependent  incompressible  viscous  Navier-Stokes  equations 
for  two-dimensional  flows  are  presented.  The  primitive  variables  of  fluid 
velocity  and  pressure  are  used.  The  Reynolds  averaged  Navier-Stokes  equa¬ 
tions  are  obtained.  Boundary  and  initial  conditions  are  given.  A  modified 
turbulence  model  based  on  Prandtl’s  mixing  length  theory  is  described  which 
expresses  the  Reynolds  stresses  in  terms  of  mean  flow  quantities  for 
separated  adverse  pressure  gradient  flows.  An  expression  for  the  force 
exerted  by  the  flow  field  on  the  body  is  derived.  The  equations  and 
boundary  conditions  are  transformed  by  the  curvilinear  transformation 
discussed  in  Section  II  and  written  in  the  computational  plane. 


A .  Basic  Equat ions  of  Motion 

The  time  dependent  incompressible  Navier-Stokes  equations,  which 
describe  conservation  of  linear  momentum,  in  orthoncrmal  indicial  noLation 
are  given  by 


-  +  s . ( V . V  . ) .  =  — c.p  +  3.  T  .  . 

it  3  j  i  i  J  Ji 


(27) 


where  .  is  the  fluid  density,  p  is  the  fluid  static  pressure,  Y  are  the 
components  of  fluid  velocity,  and  T  ^  is  the  viscous  stress  tensor.  The 
subscript  i  =  1  corresponds  to  the  x  dirction  and  i  =  2  to  the  v  direction, 
Conservation  of  mass  for  an  incompressible  flow  is  given  by 


3.  V.  =  0  (281 

3  3 


be 


Turbulent  flow  is  introduced  by  assuming  that  the  flow  quantities 
described  in  terms  of  mean  and  fluctuating  quantities  in  the  form 


V  .  =  V  .  +  V  . 

ill 


p  =  p  +  p’ 


t .  =  t . .  + 

Ji  Ji  Ji 


can 
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drop  the  viscous  divergence  tern,  and  keep  the  convective  divergence  ter::. 


as  a  correction  term.  Thus,  the  appropriate  set  of  equations  becomes 


— -  +  9.  (V.V,)  =  -  3  .  p  +  (9 . 3  .  V  .  ) 

ct  j  j  i  i  Re  '  j  j  i' 


D  =  3  .  V  . 
J  J 


(38) 

(39) 


The  calculation  of  pressure  in  the  primitive  variables  formulation 
must  be  carefully  handled.  The  pressure  does  not  appear  in  the  continuity 
equation  nor  does  a  time  derivative  of  pressure  occur.  Two  approaches 
have  been  extensively  used  to  overcome  this  difficulty.  Both  techniques 
compute  the  pressure  by  utilizing  the  velocity  divergence  I)  as  a  correction 
factor.  The  calculated  pressure  is  made  to  satisfy  continuity.  The  first 
method,  introduced  by  Chorin  (6S),  computes  pressure  using  an  iteration 


procedure  with  D  as  the  correcting  term.  The  expression  is 
(s+1)  _  (s) 


=  p 


I  p 


(in) 


where  s  denotes  the  iteration  number  and  '  is  an  acceleration  parameter 
which  can  be  either  a  constant  or  can  be  related  to  a  successive-over- 
relaxation  (SOR)  solution  of  a  Poisson  equation  for  pressure  as  shown  by 
Hodge  (53).  The  second  method  uses  a  Poisson  equation  for  pressure  derived 
from  the  divergence  of  the  momentum  equation  expressed  in  Equation  38, 

The  form  of  the  equation  is  shown  in  Equation  41. 


D  =  7  •  R  - 
t  - 


(41) 


where  Ft  represents  the  convective  and  viscous  terms.  If  the  R  expression 

is  algebraically  simplified  with  D  used  wherever  possible,  the  following 

simplified  incompressible  Poisson  pressure  equation  is  obtained 

t  a 

ft  =  -  .  <  up)  +  (vP)  +  u"  +  2v  u  +  v~  +  up  +  vP 
t  X  V  X  X  V  V  x  V 


r—  (D  +  D  )  -  (p  +  p  ) 
<e  xx  vv  Kxx  yy 


+ 


This  equation  contains  both  spatial  and  time  derivatives  of  the  velocity 
divergence  as  well  as  the  fluid  velocities  and  pressure.  Again,  continuit 
is  used  as  a  correcting  factor  for  calculating  pi  .  ssure.  Further  simp. lit i 
cation  can  occur  if  small  variations  in  the  velocity  divergence  are  as  sum*. 
In  this  case,  the  spatial  derivatives  of  1)  are  set  to  zero  in  Fquation  42 
and  the  result  becomes 

i  ■■> 

D  +  (u"  +  2v  u  +  v~)  =  -  (p  +  p  )  (4  3 

t  X  x  v  y  xx  *yy 

Hodge  (53)  used  both  Equations  42  and  43  and  found  no  apparent  difference 
in  the  flow  field  solutions  for  his  case  of  laminar  flow.  Also,  Chorin's 
method  of  Equation  40  is  related  to  the  pressure  Equation  43.  Hodge  (53) 
showed  that  the  Chorin  technique  is  approximately  related  to  a  solut  ion 
of  the  Poisson  pressure  equation  using  SOR  iteration  on  a  coarse  grid  of 
twice  the  spacing.  The  Poisson  pressure  equation  technique  is  chosen 
for  the  interior  flow  field  because  of  the  increased  coupling  which  occurs 
in  a  finite  difference  representation.  The  body  surface  pressure,  how¬ 
ever,  is  calculated  using  the  iteration  technique  of  Equation  40.  Here 
mass  conservation  is  imposed  directly  and  the  difficulty  of  evaluating 
the  Laplacian  of  pressure  at  the  hod;  surface  in  Equation  43  is  avoided. 

The  following  set  of  equations  taken  from  Equations  38,  3°  and  4* 
becomes  the  system  used  for  the  solution  of  incompressible  turbulent 
viscous  two-d imens iona 1  flow  over  airfoils: 
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numer  i  ca  1  solution  ot  the  governing  equaticn 


forcing  the  computations  in  the  computational  plane  o!  t. lined  v.  i  t  r .  t 
vilinear  transformation  discussed  in  Section  '1.  11. is  -  et  hod  reqi.ir 

the  equations  art.-  transformed  to  the  computat  ion  a  plane.  "hi  cover 
equations,  Kquatiorts  44  through  47,  art  t  rans :  ormed  u-i:.g  the  r<  la:  i 
in  Section  11  and  Appendix  A.  The  resulting  set  <  f  c  rr<-;  ..  ::d  t  r 
formed  equations  hecor.es: 


u  .  ,  v  , 

u.j.  +  —  (y,  u  r  -  y  •  u  )  +  q  ( .  u_ 


1  ,  1  ( 
j  <>'.  p  .•  1>.  >  +  kl: 


+  uf  (---)  +  u  J  ,  . 

,r  '  ,r 


VT  +  7  (yr  v  e  -  v  .  vr )  +  '■  (:<  .  v,  -  >q  v  .)  + 


~  7  (>:  :  P.  -  X.  I'  P  +  p.;; 


1  ( 


+  V,  t  1  + 

.r  ,i" 


D.J  +  •  ;  ( u  .  -  y  .  u  1 '  +  .'  (  ■  v 


+  (  v  .  V  -  X  V  .  ' 


P.  «  +  P  ■(  . 

r  .r 


=  <• 


+  (>: 


h.  Boundary  an>!  Initial  Conditions 

The  complete  definition  of  the  flow  problem  requires  the  spe,  i :  it  at 
of  sufficient  boundary  and  initial  conditions.  The  Nav ier-S t ekes  equation, 
are  coupled,  nonlinear,  first  order  in  time,  second,  ,rd,-r  in  s.pa,e,  ellipti 
partial  differential  equations.  The  elliptic  nature  of  the  equations  for 
velociiv  requires  boundary  conditions  at  each  boundarv.  1  he  first  order 
derivatives  of  pressure  require  a  boundary  condition  describ  ir.  .•  the  f  ree- 
st  rear,  pressure.  The  time  derivative  terns  require  initial  condition  on 
velo.  itv  and  pressure  throughout  the  flow  field. 


•unuarv  c; 


and  it  ions  on  the  airfoil  surface  are  con-  id, r,d  •  ir  :  . 


con  t  i  nun:-,  no  slip  boundary  condition  for  a  viscous  fluid  at  a  st.it  io;:.,r- 
solid  surtucc  is  imposed  or.  the  tangential  Velocity  at  t  h,  surface.  In 
ad..!  i  t  ion ,  the  nor". a!  component  of  velocity  must  he  zero  ‘‘or  a  stat  i.  n.ir 
surface  wit:,  no  transpiration.  These  conditions  are  met  by  tin  specif  icati 
C>:  the  !'i  rich  let  boundarc  conditions 


v,  ~  C 


w’.i  re  tin  sod.  si  ript  1)  den,  t,"  bod',  surface .  in,  pnssur,  on  t  h,  a  i  n  o  i 
surface  i-  computed,  and  in  bed."  boundary  condition  is  ft  .uirtd  by  tin  fir 
derivative  pressure  ter:  ;. 

far  fit  Id  boundary  condition--  for  velocity  and  yrisi-ur,  ,-p ,  .  ;  :  •. 
tin  frec-t  rear,  f  !  ov  :  i  ■  Id.  let  r  in  a  position  vector  in  tin  flow  f  it  !d 
wit'  tin  origin  at  the  airfoil  mid-t  ;.«*r  ! .  The:.,  a  t  h,  distance  r  !  iv:v  t! 


.nil’ll  n 


lur.:,  tiie  fluid  velo,  ity  and  present,-  approach  tin  ir 


non,!  i  mere  iona  !  f  retst  rear  value 


|  i  -t  sir.  . 


w:,<  r.  i-  r an/U  «•!  at  t«.  k  .•!  t !.«  airfoil  chord  wifi,  respt-  t  t '  1  • 

1  rei  t  r>  a:.. . 

The  f  ir<t  orji  r  t  ir..c  derivative  terms  require  initial  eond  i  t  ion-  r 
tin  vi'i  ,  it v  an.!  prcs.-ur,  in  the  f  ield.  Fro:  a  t  ii  a1,  viewp.  i:.t  , 

since  tv.  pt  tic  .!,  jf  s.-’.ut  ion  i  s  !es>  i  red  ,  ;»Tiv  .  or!  inu.-  J  •.  i  r:  :  t  mil  s;  .  .  :  ’  . 

eat  i  ■  o:  Vi  '..a  it  -  an.!  pressure  (  subject  to  the  prescribed  '  ■  •  i : .  o  a  r ;  end  i  t  ’ 

i  s  p  c  rm  i  s  s  i  1  ! .  . 

I  in  bound,  ir'  c  ond  i  t  ions  and  initial  rent!  it  ions  an  al  s„  t  r.uis:  ern.e  J  to 
the  comput at  ionu!  pi  ant  .  ;  !.<•  values  lor  ear!,  conditi-  •  remu  ir  unal  te-red  art! 

an  thr-  :  repeated  !:t  re  .  !  :.«•  lot  at  ior.  in  : t  ran- :  armed  pi. on  depend-  .  r. 

the  s ; 1 1* .  i:i.  t  run.-:  err  .at  ior.  >  ho..:,  to  generate  t!.e  comp  utat  ion...’,  grid.  Th  i 

ratter  i  J  iscus.se.!  in  Section  IV  where  the  numerical  procedure?  on  descr  i’-cd  . 
'.he  setts  it  iv  it  v  analysis  of  the  outer  boundary  location  is  addressed  in 
Appendix 

f.  Turbulence  V.hK_1 

Turl-ult  r.ce  model  1  ing  is  of  considerable  interest  in  theoretical,  experi¬ 
mental  ,  and  con.pulat  iona  1  research,  wort.  Matty  .-u  rodynamic  problett.s  of 
pra.  t  ical  interest  contain  regions  of  turbulent  flow.  Several  model,  have 
been,  developed  which  vary  in  complexity  from  Prandt  1  ’  s  mix  in.-  let  thief; 

to  techniques  employing  auxiliary  differential  equation  for  mix!::/  len/t: 
viscosity  coefficients,  and  turbulent  kinetic  e:n  r/ies.  The  a.vur.e  ;  of 
each  model  varies  with  both  the  geometry  and  the  flow  i  i.  Id  quart  itv  .:  let 

comparing  results.  In  computational  work,  tht  abi  1  it-  t  o  t  ra:..  !  a! .  tur¬ 

bulence-  model  into  numerical  procedures  and  the  computer  re-  ur.c  r»  .u  i  i ,  ,! 
are  additional  factors  which  need  to  he  cons i der i d . 

In  this  investigation,  turbulence  is  node  !lu!  with  an  a  ! 1  r..  . .  id.: 
viscosity  approach  based  on  Prandt  1  *  s  mix  in/  length  t!.o>r.  .  1 :« 

concept  of  this  method  is  that  the  turbulent  ke-.nold  ;,t  r.  >«  ■  ,  -  V  .  . 


car.  he  related  to  the  strain  rate  of  the  mean  flow  in  a  wav  analagous  to 
that  for  laminar  viscous  stresses.  The  fundamental  expression  was  previous!;, 
introduced  in  Equation  3d  to  obtain  the  simplified  turbulent  incompressible 
Nav ier-Stokes  equation,  Equation  33.  The  algebraic  eddv  viscositv  has 
been  used  extensivelv  in  obtaining  numerical  solutions  for  a  variet'  of 
aerodynamic  problems  (24,25,p9) . 

The  eddy  viscosity  model  for  turbulent  flow  with  zero  pressure  gradient 
near  a  solid  surface,  which  incorporates  the  mod i f icat ions  of  Cebeci  and 
Smith  it-il.  Van  Priest  (13),  and  Peiwert  (57),  is  given  below  and  then  discus 


(k,  vP,  1  (i r  +  v“) 

11.  V  V 


=  1  -  exp 


(---—) 


1  +  5.5  (v/:) 


(1  -  7-I  dv 
u 


The  tangential  surface  direction  is  given  by  x  and  the  normal  surface 

direction  bv  v.  The  inner  and  outer  eddv  viscosities,  and  ,  mod, 

i  o 

tin  I  aw  of  t  he  Wall  and  law  of  the  Wake  (70)  regions.  The  inner  rr.ixir.  ,■ 
length  is  assumed  proportional  to  the  distance  v  from  the  surface.  k. 
the  Von  Kurman  constant.  The  Van  Priest  (.I’D  damping  factor  nodi  Is 
the  laminar  sul-luver  region  where  •  is  the  surface  shear  stress  and  i* 
the  decay  constant  .  The  d i -placement  thickness  '  de:  ined  bv  !‘<:uat  ten  : 
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t  oni 
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nondiriensional  iz<?d  by  tlu*  local  boundary  layer  distance  .  The  be.;! 

edre  boundary  lavor  velocity  u  is  the  outer  l aver  characteristic  veloc 

-  t. 

1  he  standard  Cebec i-Smith  (63  >  values  lor  the  constants  are  kj  =  .41 , 
k,  =  .0!  >  and  A  =  2b.  Continuity  of  the  eddy  viscosity  distribution, 
achieved  by  sv  it  chine  f  ron.  the  inner  to  the  outer  n.odel  whenever  .  fir 
exceed. 

Recently,  do’e  and  Hankey  (71)  compared  this  eddy  viscosity  mode! 
bv  Kquat  ions  55  through  5°  with  experimental  data  for  several  constant 
adverse  pressure  gradient  attached  flows.  Thev  varied  the  k, ,  1.,,  and 
parameters  in  a  numeric  si  turbulent  boundary  lay*,  r  program  until  t  he  s, 
tion.s  matched  the  experimental  boundary  layer  velocity  profiles  and  ski 
friction  ooe  f  f  ic  ic-r.t  s  .  These  parameters  had  substantially  different  va 
than  the  previously  cited  values  obtained  for  zero  pressure  gradient  : 1 
They  found  that  when  the  flow  encounters  an  adverse  pressure  gradient  t 
eddy  vinco.-ity  parameters  rapidly  adjust  to  the  values  of  k,  =  .65  and 
A  =  52.  The  outer  eddy  viscosity  constant  appears  to  linearly  increas 
wit:,  distance  measured  from,  the  start  of  the  adverse  pressure  gradient  . 
Thus  .lobe  and  Hankey  (71)  found  that  k  ,  car.  vary  from:  a  value  of  .  P 1  6 c 
value  s  of  .03  or  greater  given  that  the  adverse  pressure  gradient  exist 
for  a  sufficient  length  downstream. 

The  use  of  an  eddy  viscosity  model  similar  in  form  to  the  one  give 
in  Equnt  ions  55  through.  59  for  an  airfoil  at  angle  of  attack  require.- 
examination  of  the  typical  flow  field  characteristics.  The-  upper  surf  u 
of  the  airfoil  has  a  pressure  distribution  which  includes  a  region  of  1 
pressure,  tin-  suction  peak,  near  the  airfoil  leading  edge.  The  sik  t  i  • 
is  followed  bv  an  adv<  r.  i  pressure  gradient  region.  The  pres.-uri  ,  r  .  1 

tends  to  decrease  a  the  flow  recovers  toward  the  trailing  edt.  ■ 


11 


turbulence  nod el  that  is  used  for  an  upper  surface  turbulent  ban. car 
should  include  the  effects  of  this  adverse  pressure  gradient.  "be  in 
eddy  viscosity  parameters ,  therefore,  take  on  the  Jobe  and  Han he v  'll 
of  k1  -  -bS  and  A+  =  52  in  this  region.  The  outer  edd"  viscu.-  if.  par 
k.-,  mav  also  change  in  the  adverse  pressure  gradient  region .  However , 
rate  of  change  is  affected  bv  the  pressure  gradient  and  downs! rear,  di 
The  pressure  gradient  does  decrease  rapidly  as  the  flow  proceeds  down 
Any  tendency  for  k ,  to  increase  is  probably  offset  by  the  decreasing 
pressure  grad  tent .  Thus,  the  value  of  k ,  is  kept  at  the  ecutli;  riur. 
of  .  (!1  (-;s  . 

The  eddy  vise. -si  tv  mod.  ’.  with  the  discussed  adverse  pressure  gra 
mod  if  icat  ions  is  now  written  in  ter  r..-  of  the  nondimens  ional  variables 
duced  in  liquations  S-*  and  37.  All  length  quant  itios  arc  nondir.tnsior. 
with  respect  to  the  airfoil  chord.  The  mean  and  nondimons ional  quart 
are  henceforth  implied.  The  eddy  viscosity  model  becomes 


1  ■) 


=  (k.vh  (u"  +  v")  Kc 
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(1  -  1  dv 
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when 


+ 


<■:  ar.c.  a  are  considered  adverse  pressure  gradient  para". etc 

ck’v  viscosity  node]  must  also  represent  the  phvsical  eharac 


ii.  larvc  separated  regions 


Tra i 1  in 


which  occur  during  airfoil  stall. 


stall  occurs  when  the  turbulent  boundary  layer  separation  roves  toward  the 
leading  edge  of  the  airfoil  as  the  angle  of  attack  increases.  The  pressure 
distribution  in  a  trailing  edge  separated  region  has  snail  gradients  (72). 

This  observation  motivates  another  change  to  the  eddy  viscosity  model.  The 
inner  eddy  viscosity  parameter  is  relaxed  back  toward  the  equilibrium  value 
of  .41.  The  separation  point,  which  precedes  the  return  to  zero  pressure  grad¬ 
ient,  is  selected  as  the  point  to  start  the  exponential  decay.  The  following 
relaxation  form  with  distances  nondinensionalized  by  the  chord  is  used 

s  -  s  , 

f  =  A  1  +  B  exp(-  — — )  s  -  s 

s  d 

r 

('■:■>> 

f  =  1  s  s  . 

d 

n 

where  f  is  the  relaxation  factor  which  multiplies  the  initial  value  of  k“, 

s  is  the  distance  downstream  from  the  trailing  edge  separation  point, 

s,  is  the  delav  distance,  s  is  the  relaxation  distance,  and  A  and  B  are 
d  ‘  r 

constants.  The  constants  A  and  B  are  determined  by  applying  the  conditions 

2  0 

f  =  1  at  s  =  and  f  =  k^/k^  as  s  approaches  infinity  where  is  the 
asymptotic  value  of  k-j  downstream.  Thus,  A  and  B  hecone- 

A  =  k“f/k^  B  =  A_1  -  1  (*6) 

Measurements  of  Bachalo  (73)  and  Baker  (74)  in  the  trailing  edge  separa¬ 
tion  region  indicate  the  presence  of  a  lower  but  finite  turbulent  intensitv 
when  compared  with  the  separated  turbulent  shear  layer.  The  velocity  deriv¬ 
ative  terms  in  the  inner  eddy  viscosity  and  damping  factor  expressions,  Equa¬ 
tions  60  and  62,  can  approach  zero  in  a  region  near  a  velocity  inflection 
point.  In  this  case,  a  zero  eddy  viscosity  is  permitted  in  a  turbulent 
separated  region  which  is  contrary  to  the  apparent  flow  structure.  This  same- 
zero  eddy  viscosity  behavior  can  als-o  occur  in  the  case  of  sheddine  separation 
bubbles  from  the  airfoil  surface.  In  order  to  prevent  the  appearance  of 
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laminar  regions  within  the  turbulent,  adverse  pressure  gradient,  separated 
flow  field  during  airfoil  stall,  a  limiting  technique  is  employed.  The  inner 
eddy  viscosity  is  prevented  from  decreasing  in  value  as  the  normal  distance 
from  the  surface  increases.  This  restriction  simulates  the  uniform  turbulent 
intensities  measured  (73,74)  in  the  trailing  edge  separated  region.  The  innc r 
eddy  viscosity  is  further  limited  by  imposing  a  condition  of  no  decrease 
in  the  downstream  direction.  This  limit  provides  for  a  finite  turbulent 
intensity  in  separated  regions  which  may  develop  in  the  turbulent  boundary 
layer.  The  limit  is  initiated  by  obtaining  the  distribution  of  eddy  vis¬ 
cosity  in  the  attached  boundary  layer  near  the  leading  edge  of  the  airfoil. 

This  distribution  C  is  compared  with  the  distributions  calculated  by  the 
model  at  downstream  locations.  The  outer  eddy  viscosity,  which  has  no  velocity 
derivative,  is  not  limited  in  any  way.  The  "limiting"  technique  is  passive  in 
the  sense  that  the  locally  computed  value  for  the  inner  eddy  viscosity  is  used 
whenever  possible.  Numerical  experiments  without  the  "limiting”  technique 
displayed  unphysical  results  because  the  conventionally  modelled  eddy  viscosity 
is  divergent.  The  inception  of  separated  flow  decreased  the  eddy  viscosity 
thereby  causing  futher  separation.  The  "limiting"  concept  is  necessary  to 
prevent  this  unphysical  result. 

The  transition  from  laminar  flow  to  fully  turbulent  flow  is  modelled 
using  the  transition  model  of  Dhawan  and  Narasimha  (75).  The  expression  for 
the  transition  factor  r  is  given  below 
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model  are  conveniently  expressed  in  terms  of  quantities  in  the  transformed 
plane  described  in  Section  II.  The  velocity  derivatives  in  the  inner  eddy 
viscosity  model  are  formulated  by  determining  the  tangential  and  normal 
components  of  velocity,  u  and  v  respectively,  to  the  body  surface  r  =  ’  j 
at  a  position  ’ 

u(  f,  ’■)  =  c9u(  r,  r)  -  c^v(f,  ’  )  (bb) 

v(  r,  h)  =  c1u(  i,  r)  +  c,v(  f,  r;)  (69) 

where  and  c0  are  components  of  the  unit  outward  normal  to  the  body  sur¬ 
face  at  position  (  F,  r ,.)  given  hv 


c.  =  -v  ./>  •>  c  =  — 

1  f  -  ;  >  t 

Then,  the  directional  derivative  definition  is  used  to  obtain  the  surface 
normal  derivative  of  u  and  the  surface  tangential  derivative  of  v  expressed  as 


-  u  t 

-  Ci(>iU 


-VOJ 


)  +  c9  (:■;  -  x,.u  .)  /  J 


(70) 


~  =  c,(v  v  -  y  fv  )  -  c  (>;  o'  -  x  v  .) .  /J  (71) 

■  s  '  T,  i  >  ri  it’  r,  ■ 

which  correspond  to  the  u  and  v  derivatives  of  the  Cartesian  eddv  vis- 
r  y  x 

cosity  model  expressions.  Equations  60  through  64. 

The  turbulence  structure  in  the  far  wake  is  modelled  with  an  eddy 
viscosity  expression  similar  to  the  outer  eddy  viscosity  model  in  the  tur¬ 
bulent  boundary  layer.  The  model  is  based  on  the  assumption  of  a  self¬ 
preserving,  equilibrium  turbulent  wake  with  co-stant  pressure  in  incompres¬ 
sible  flow.  The  approach  has  been  experimentally  investigated  by  Narasimha 
and  Prabhu  (76)  who  obtained  the  following  expression 


=  .  k.,w 
w  i  o  w 


(721 


where  w  (x)  =  max  ■ b  -  u(x,v)  is  the  maximum  velocity  defect  m  the  wake, 
o 

v 

1  is  a  half  wake  width  defined  where  the  velocity  defect  becomes  one  half 
w 

the  maximum  value,  and  k,^  is  a  constant  equal  to  .065.  If  the  t-xper  inter,  la  1 1 ; 
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integrated  value  of  the  wake  defect  function  is  used,  then  Equation  72  can 
be  written  in  terms  of  the  parameters  of  edge  velocity  and  displacement 
thickness  similar  to  boundary  layers  as  follows: 

—  =  k,u  6*  Re  (73) 

L.  3  e  v 

where  u  is  the  nondimensional  velocity  at  the  edge  of  the  wake,  is  the 
e  w 

nondimensional  displacement  thickness  for  the  half  wake  as  defined  by  Equa¬ 
tion  64,  and  k^  becomes  .0634.  This  expression  has  been  utilized  by  Green, 
et  al  (77)  in  an  integral  method  for  predicting  two-dimensional  incompressible 
and  compressible  turbulent  wakes  of  lifting  airfoils.  For  an  incompressible, 
constant  pressure  wake  the  momentum  integral  equation  reduces  to  a  statement 
that  the  momentum  thickness  is  a  constant.  These  conditions  exist  in  t 
far  wake  of  an  airfoil.  In  addition,  in  the  far  wake  the  displacement  d 
momentum  thicknesses  become  equal  as  shown  by  the  experimental  r.  ,ults  of 
Narasimha  and  Prabhu  (76).  Thus  the  far  wake  eddy  viscosity  is  approximate)'- 
constant.  The  transition  in  the  wake  from  the  turb-ilent  core  region  1.0  the 
inviscid  freestream  region  on  each  side  of  the  wake  is  modelled  with  the 
intermittent;.-  factor  of  Equation  63  where  f  is  now  the  wake  half  width  and 
y  is  the  normal  distance  from  the  wake  center. 

The  turbulent  near  wake  region  mixing  rates  where  the  wall  boundary 
layers  merge  and  eventually  become  the  far  wake  are  not  well  understood  fro:;: 
either  an  analytical  or  experimental  viewpoint.  In  the  flow  regime  under 
consideration,  a  laminar  boundary  layer  from  the  lower  airfoil  surface  mixes 
with  a  turbulent  boundary  layer  from  the  upper  surface.  An  edd\  viscosity 
model  similar  to  that  used  in  the  outer  turbulent  boundary  layer  a:u!  the  far 
wake  is  used.  Inouye,  et  al  (78)  applied  the  Cebeci-Fmith  (A3)  outer  layer 
model  in  the  near  wake  of  a  flat  plate  and  obtained  good  agreement  for  the 
mean  velocity  profiles.  They  allowed  the  constant  k  to  reach  twice  t 


equilibrium  value  of  .0168  over  a  length  of  one  chord  from  the  trailing 
edge.  The  computed  eddy  viscosity  was  still  about  50'  low  when  compared 
with  experiment  at  the  one  chord  distance.  Later,  Burggraf  (79)  compared 
several  turbulence  models  for  various  near  wake  flows  and  concluded  that 
the  Cebeci-Smith  model  adequately  predicted  the  mean  velocity  prefiles. 

The  merging  of  the  upper  and  lower  surface  boundary  layers  is  modelled 
with  a  simplified  interaction  hypothesis  approach  introduced  by  Bradshaw  (80). 
The  two  layers  are  initially  allowed  to  develop  downstream  a  distance  WL^ 
without  interaction.  At  this  point  the  eddy  viscosity  is  increased  exponen¬ 
tially  in  the  adjacent  laminar  boundary  layer  from  a  zero  value  to  the  far 

wake  calculated  value  c  at  a  distance  WL„  from  the  trailing  eu;e.  This 

ow  2 

eddy  viscosity  variation  is  formulated  as  follows 

s  -  WL 


c  =  1/2  e  . 1  +  tanh 
ow 

•:  =  0 


(8 


WL0  -  WL1 


): 


WL, 


s  '  VL1 


(79) 


where  WL  t  is  the  distance  from  the  trailing  edge  to  a  loeatien  half  way 
between  the  distances  WL^  and  WL0  and  s  is  the  local  distance  from  the 
trailing  edge.  The  coefficient  8  is  selected  so  that  •  0  at  s  =  WL^ 

and  ■:  -  c  at  s  =  WL.,  because  at  s  =  WL0  the  hyperbolic  tangent  function 
becomes  equal  to  0.9993  while  at  s  =  WL ^  it  becomes  -  0.9993.  The  eddy 
viscosity  profile  in  the  turbulent  boundary  layer  near  the  trailing  edge 
is  extended  into  the  near  wake. 


D .  Force  Coefficients 

The  force  which  is  exerted  on  the  body  by  the  moving  fluid  around  the 
body  can  be  determined  from  the  Cauchy  integral  equation  for  conservation 
of  linear  momentum  applied  at  the  body  surface  or  on  a  closed  path  surround¬ 
ing  the  body.  The  basic  expressions  for  the  force  coefficients  formulated 
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•~2x  rp 


where  the  vorticity  on  the  body  surface  is  given  by 


-  (v  .v  +  X  Ai  )  /  J 


The  corresponding  moment  (positive  counterclockwise)  about  a  point  P 

located  within  the  bodv  cross-section  with  coordinates  (x  ,  v  )  is  given 

P  '  P 
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(x  -  Xp)(-2xrp  -  ^— )  -  (y  -  yp)  (-y  rP  -  ~j~ 


-)  d  (H2I 


Next  consider  again  the  more  general  case  of  a  closed  path  around 
the  body  whose  force  coefficients  are  given  by  Equations  P.12a  and  D.ldb. 
Assume  that  the  closed  path  is  an  r  contour.  Then,  using  the  transformed 
expressions  in  Appendix  A  for  the  outward  normal  vector  to  an  ’  contour, 
integrals,  and  derivatives,  the  force  coefficients  for  the  body  in  term.-. 


of  quantities  on  a  constant  ’  line  in  the  flow  field  become 
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where  r),  <•  r;  <  q  .  The  terms  in  the  line  integral  represent  in  order 
b  —  p  —  max  r 

pressure  forces,  viscous  forces,  and  convective  outflows  of  linear  momentum. 
The  area  integral  term  is  the  time  rate  of  increase  of  linear  momentum  con¬ 
tained  in  the  control  volume  bounded  bv  the  bodv  T,.  contour  and  the  selected 

n 

ti  contour. 

P 

The  corresponding  lift  and  drag  coefficients  are  then  calculated  using 
Equations  D.13  and  P.14. 


40 


SECTION  IV 

NUMERICAL  METHOD 

T)it-  numerical  methods  used  to  obtain  solutions  for  im  «c.j  res-  il  e 
turbulent  flow  over  airfoils  are  presented.  The  numerical  f  ini  It  d i ; tere 
procedure  which  yields  the  grid  transformation  is  first  discussed.  The 
implicit  finite  difference  method  for  obtaining  approximate  solutions,  tc 
the  unsteadv,  incompressible,  turbulent  Navier-Stokes  equations:  is  next 
described.  Numerical  boundary  conditions  and  initial  conditions  are  con¬ 
sidered.  The  numerical  procedure  for  incorporating  the  turbulent  eddy 
viscosity  model  is  discussed.  The  computation  of  the  force  coefficient- 
on  the  bod;.'  is  then  presented. 

A  .  G  r  i_d  J  r  a  n sfo  rmat  ion 

The  numerical  solution  of  the  governing  equations  for  t  he  t  nir.sf  or::  a 
tion  is  carried  out  in  the  transformed  or  computational  plane  on  a  square 
mesh  (  -)  s vs tern.  The  bodv  contour  C,  and  the  far  field  houndarv  < 

become  constant  ’  contours  (Figure  1).  The  far  field  boundary  is  usually 
used  to  approximate  infinity  in  incompressible  flow  problems ,  Constant  1 
contours  in  the  physical  plane  are  required  to  transform  to  equi-s'pa.'od  ' 
coordinate  lines  with  unit  step  size  in  the  transformed  plane.  The  spa, i 
contours  in  the  physical  plane  is  controlled  by  the  set  of  elliptic  equal 
The  ■  contour  spacing  on  the  body  surface  and  outer  boundary  in  the  jd.vsi 
plane  is  determined  by  the  chosen  point  distribution  on  each  surf  act  .  1 h 

constant  •  contours  are  also  required  to  transform  to  equ i -spaced  "  coor¬ 
dinate  1 ines  with  unit  step  size  in  the  transformed  plane.  The  line: 
are  denoted  by  the  subscript  i  with  range  1  •_  i  '  ] MAX  where  IM.V-  is  the 

of  ■  lines.  Similar iv,  the  ’  lires  are  identified  by  a  subscript  j  with 


i  on 
,  ;il 


iiu 

run 
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of  1  ■  j  •  JMAX  whore  JMAM  is  t  ho'  number  ol"  ’  lines.  '!  ho  value;  j  =  1 
and  j  =  .IMAM  art-  the  bodv  and  outer  boundary  surfaces,  respect  ivc  1  \  . 

Tlio  constant  •  lines  i  =  1  and  i  =  1  MAM  denote  the  branch  cut  in  the 
physical  plane  and  are  therefore  identical. 

The  governing  elliptic  transformation  equations,  liquations  ua  and  °h , 
are  written  with  finite  differences  in  quasi-1 inear ized  fort.  fur  the  loc¬ 
ation  (i,  j).  Second  order  accurate  central  differences  given  in  Append::-: 
are  used  for  evaluating  each  term.  The  highest  order  derivative  contain, 
the  unknown  at  location  (i,  j)  while  the  transformation  coefficients,  whic 
contain  lower  order  derivatives,  are  evaluated  with  previous  values  of  t 
variables  x  and  y.  The  difference  expressions  are  piv«.r.  below  where  the 
indices  i  and  j  are  understood  when  omitted. 


x  .  .  =  a  ( x  .  ,  .  +  x  .  .  )  -  sr  ( x  .  .  .  ,  ,  -  x .  ,  ,  .  ,  +  x  .  ,  .  ,  -  x  .  ,  .  ,  .  1 
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:(yl+l  +  yi-l)  2  (yi+l  j+1  y i+ 1  j-1  +  y i - 1  j-1  ~  vi-l  j+15 


+  +  +  r  f(vi+i  •  yiVF  +  fyi+i  -  yi-P°:  /-(- +  •’ 


The  i  index  near  the  branch  cut  is  adjusted  as  follows: 

when  i  =  1  then  i  -  1  -  IMAM  -  1 

when  i  =  IMAM  -  1  then  i  +  1  =  1 

when  i  =  IMAM  then  i  =  1 

The  set  of  finite  difference  equations  has  the  range  of  indices  given  by 
1  -  i  •'  IMAM  -  1  and  2  '  j  p  JMAM  -1  which  results  in  ( 1  MAM- 1 )  (.IMAM- 2 1 


equations  for  the  unknown  values  of  x  and  y  at  the  interior  grid  points. 
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I  he  system  of  difference  equations  is  solved  using  successive  over  relaxa¬ 
tion  (SCR)  iteration.  For  a  general  function  f,  an  approximate  solution 
for  f  using  SOR  iteration  has  the  form 


f(srl)  =  f~  +  (1  _ 


;  i-  / ) 


where  s  represents  the  iterate  number,  is  an  acceleration  parameter,  and 

*  represents  the  current  value  of  f  calculated  from  the  difference  equation. 

Tile  acceleration  parameter  is  always  one  for  the  first  iteration.  The  system 

of  linearized  difference  equations  given  by  liquations  85  and  86  obtained 

frot.  the  elliptic  generating  differential  equations  is  consistently  ordered  ( - 

ihus  local  optimum  acceleration  parameters  ...  can  be  calculated  fro::,  .■-.air  i 

tj 

theory  (81)  and  are  given  by  Hodge  ( 5?) .  Convergence  of  the  iteration  procedure 
is  determined  by  comparing  the  absolute  value  of  the  difference  in  successive 
iterates  with  a  prescribed  error  criteria  V. 

The  SOR  iteration  technique  requires  an  initial  guess  for  the  solution. 
Geometric  contours  with  a  shape  similar  to  the  outer  computational  boundary 


are  used,  for  constant  ’  lines  and  straight  lines  which  emanate  frot.  tin  bod; 
center  are  used  for  ’  lines.  A  similar  approach  has  been  successfully  used 
by  Thames  (44)  and  Hodge  (631. 

The  forcing  function  0(’,  ’1  defined  by  liquation  13  is  evaluated 
throughout  the  (  \)  plane  following  the  procedure  described  in  Section  1  ]  .  P. 

The  similaritv  parameter  values  '  (which  give  equal  increments  of  tangential 
velocity  for  a  specified  number  of  boundary  layer  points)  are  calculated  fro: 
Equation  18  using  Newt on-Raphson  iteration.  The  equivalent  v  value-  a 
given  >:  (  line)  art  found  from  the  similarity  relation,  Equation  in. 
r;  derivatives  given  l»v  Equations  22  and  2'i  arc-  then  approximated  using  the 
above  values  for  v  and  the  specified  ’  differences  as  follows: 
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where  =  1,  y.  are  the  boundary  layer  y  values,  and  '  y  is  tin  spec  i  f  i  *  d 
far  field  v  difference  between  the  last  two  points  on  the  giver:  '  line- . 

Then  tile  damping  factor  D  is  found  using  liquation  24 .  Next  tin  analytical 
expression  for  d'  /dv  in  Equation  21  is  evaluated  on  the  far  field  interval 
.‘.v  and  on  each  boundary  laver  .‘.v  interval  determined  bv  the  v.  valla  .  . 

...  •  i 

This  system  of  equations  for  A^  (  0  is  solved  using  Gauss  el  iir.ir.at  ior. .  Tims . 
for  a  given  line,  values  for  A  and  P  arc  calculated  for  use  in  eva'.uat  ir.g 
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interpolating  technique  given  in  Appendix  C . 

B .  Navi er-Sto k e V  i_n_i_u  differ e_ncc-_  Equations 

Approximate  solutions  of  the  unsteady,  Reynolds  averaged,  incompres¬ 
sible,  turbulent,  Nav ier-Stokes  equations  are  obtained  by  using  an  implicit 
finite  difference  procedure.  Explicit  numerical  methods  can.  K  u.-*  !  t  ,  • 
obtain  approximate  solutions  by  advancing  the  time  by  small  increment!  . 

The  allowable  time  increments  are  constrained  by  numerical  stability  or  it.  via. 
Imp  1  if  it  finite  difference  methods  can  also  provide  approximate  so  1  -  it  in: 
marching  in  time  with  small  time  steps.  Implicit  methods,  usual’.-  do  n  :  ...n 
numerical  stability  constraints  imposed  on  the  time  step.  iiewc  v>  r  ,  at  *  ...  1 
time  step  a  lar.n  system  of  simultaneous  difference  ccpiat  ions  rot  ’<  ■  v.  :  . 

A  matrix  it*  rat  ive  technique  such,  as  SDK  iteration  is  a::  attractive  -..Inti  - 
method  be*  ansi-  of  the  largt  matrix  si/e.  Convergence  er  i  t  <  r  i  a  ar*  1 1  :  t  t  - 

dined  in  place  oi  stability  criteria. 


The  inpl  icit  "njtluui  consists  of  hot!;  a  linearization  of  the  equat  ion.; 
ami  the  application  of  finite  differences  to  those  equations.  The  x  com¬ 
ponent,  transformed ,  momentum  equation  given  by  liquation  48  is  considered 
first.  Rc-ur rang inc  the  convective  terms,  the  equation  becomes 


I  l  i  • -  -  v 


1  -  u  (v  -  -  u  --) 


:  •  .1 


ll!  .  .  -  d .  U  .  +  .  Ut  _ 

+  v  - -  c  ■  -  +  u  -  (-— )  +  u,  (  -,i  (89) 
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liquation  S9  is  new  written,  in  quasi-linenr  form  vhe  re  lover  order 
der ivnt Ives  of  the  variables  occur  in  "coef f icients" .  The  current  v-lue 
for  the  unknown  variable  at  local  inn.  (i,  j)  occurs  in  the  highest  ore,  r 
derivative  of  each  t  r:r.  while  th.e  "coef ficients"  are  evaluated  with  tin-  lu.-t 
available  values  for  the  variables  u,  v,  and  p.  In  tb.is  form,  tin'  coefficient 
of  u.  in  tin’  convective  terms,  (uy,  -vx,  1  /.l ,  is  proportional  to  the  component 
of  Velocity  in  the  ’  direction.  Similarly,  the  coefficient  of  u.  in  the  convec¬ 
tive  terms,  (vx  .-in.  d>  /.i ,  is  proport  eonal  to  the  component  of  veloc  itv  in  the 
■  direction.  These  quantities  an  designated  IT  and  VC,  respectively. 

Finite  difference  expressions  ; iven  in  Appendix  B  are  used  to  approxi¬ 
mate  t  he  terms  in  Equation  Sy .  first  order  accurate  backward  dif f-ren.  es 
ar«.  used  for  the  time  derivative.  Tile  first  derivatives  of  pressure  are 
up;- r  i  x  im.i  t  ed  with  second  order  accurate  central  differences.  likewise, 
soc-th!  ei  r  de  r  central  difference  expressions  are  used  for  the  so,  end  dor  i  vat  i  ve  : 
or  velocity  in  the  viscous  terms.  Second  order  accurate-  upwind  differencin' 
i-  u  ed  for  the  velocif.  first  de-r  ivatlves  in  the  conveetivi  t  err..- .  1 ; ., 

quant  it  ios  "('  and  VC  determine  tin  sign  of  the  local  velocity  in  t he-  and  • 
direct  ions.  It  IT  then  backward  upwind  d  i  fft-rmring  is  used  lor  t  lit 

derivative  while  ’T  h  dictates  forward  upwind  differencing.  .-Ms,-. 
vird  -ipvfr  !  ,!i‘!i-r.'-wi:u’  for  -  derivative-  occurs  if  Yt  !>  v?  i  1 1  f  or.-. r,.! 


upwind  differencing  i.s  used  if  VC  •  0.  Second  order  accurate  upwind  d  i 

ferences  are  also  used  for  the  velocity  first  derivatives  in  the  viscous 
terms.  The  direction  is  deterr.im  d  by  the  sign  of  the  coefficient-  . 
Points  on  the  .1  =  2  and  J  =  JMAX  -  1  contours  use  fir.-t  order  upwind 
differences  for  ’  first  derivatives  if  the  velocity  direction  indicates 
flow  from,  the  body  surface  (J=ll  or  tin.-  outer  boundarv  (  J  =  i.  ;;Lar 

the  body  surface  the  grid  spacing  is  extremely  small,  and  velocitv  gradie 
are  small  near  the  outer  boundarv.  Thus,  in  both  cases  the  associated 
artiiicial  viscosity  has  a  negligible  effect.  The  transformation  deriva¬ 
tives  in  Equation  S9  are  evaluated  with  second  order  central  differer.o  •  . 
The  velocity  derivatives  in  the  velocity  divergence  term  are  also  evaluat 
using  central  differences. 

Mtferen.ee  equations  are  written  with  the  following  convenient  con¬ 
vent  ions .  Terms  which  contain  unknowns  at  location  (i,j)  are  full'.-  sub- 
scripted.  The  finite  difference  expression  for  location  (i,  ; 1  as 
described  previously  is  understood  when  only  the  term  is  given.  The  sub  - 
script-  (i,  j)  and  superscript  n  are  assumed  when  omitted.  1  he  result  in, ■ 
u  i  t  t  e  r  e  n  c  e  equation  for  Equation  becomes 


i  -  n  ~  ^ 

!!ii . +  Ciuj  -  4ulcl  +  uu"'N|  rc  +  (-5ui  ■  4u Tfi  +  u 


J  (pi+l  •  pi-l}  +  2J  (pj+I  ‘  pi-]' 


t  .J' 


T  (u  ,  -  2u .  .  +  u.  .)  +  (u  -  2i: .  . 

1+'  ij  i-l  1  1  t 


d  r 


i+l  j+l  u i-l  j  +  !  U i-l  i-l  u  i  -  !  ■-i'1 


(  ^  ’-1  -  —  -+  l 1 

1 


1 VI  +  11 1 V2 }  "V-  "  0,1  j  '  4"y  y  i 


+  i* 


f- 


where 


and  YY 
where 


for  UC  > 

•  0 

for 

rc  "  o 

IC1  =  i 

-  1 

IC1 

•=  i  +  1 

IC2  =  i 

_  7 

IC2 

=  i  +  2 

for  YC 

0 

for 

YC  '  0 

h-* 

II 

<—». 

-  1 

JCl 

=  j  +  1 

JC2  =  j 

_  7 

JC2 

=  j  +  2 

- /2  J  2 

and 

YY  : 

~r/2.T 

for  L'Y 

0 

for 

YY  •'  0 

IV 1  =  i 

-  1 

IY1 

=  i  +  1 

IY2  =  i 

_  7 

IY2 

=  i  +  2 

for  YV  > 

.  0 

for 

YY  <  0 

JY1  =  j 

-  1 

JV1 

=  j  +  1 

JY2  =  j 

_  2 

JV2 

=  j  +  2 

The  difference  equation  for  conservation  of  mass.  Equation  51,  for  location 
( i  ,  j )  is  given  by 


D  =  37  (ui+i 


U  . 
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2 J  (vi+l  Vi-1)  2J  (uj+l  Uj-1) 
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The  following  finite  difference  definitions  are  introduced  for  the 
transformation  derivatives  and  coefficients  at  location  (i,  j): 


XF.TA  =  x  /2J 

r. 

XXI  = 

X../2J 

YETA  =  y  /  2.J 
r; 

YX1  = 

y./2.T 

? 

7 

ALFA  =  '<./  J“ 

GAMA  =  y/.r 

BETA  = 

-P/2  .r 

With  these  definitions,  Equation  90  can  he  written  to  provide  an  expression 
for  un .  given  by 
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where  *  denotes  the  value  of  the  unknown  calculated  directly  from  tin:  dif¬ 
ference  equation. 

Several  observations  are  made  concerning  the  difference  equation  for 

u„  given  by  Equation  92.  The  unknowns  v_  and  p^  do  not  appear  explicitly 

In  fact  p  does  not  appear  at  all.  The  v  component  of  velocity  y  does 

appear  implicitly  in  the  quasi-linear  coefficients  1C  and  VC  as  does  u .  .  . 

^  J 

The  previous  value  for  u^.  occurs  explicitly  in  the  term  containing  the 

velocity  divergence  P.  The  finite  difference  linearization  assumption 

requires  that  the  latest  prior  values  for  u..  ami  v. .  are  used  in  these 

ij  !J 

instances . 

In  an  analogous  manner,  the  y  component,  transformed  momentum  equation. 

Equation  4°,  is  written  for  vn .  and  becomes 

tj 
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Again,  p„  does  not  appear  either  explicitly  or  implicitly.  Trie  unknown 

u_  appears  implicitly  in  the  coefficients  VC  and  VC.  Also,  v„  occurs 

implicitly  in  VC  and  VC  and  explicitly  in  the  velocity  divergence  term. 

The  most  current  previous  values  for  u..  and  v..  are  used  in  these  lower 

ij  ij 

order  coefficients. 

The  static  pressure  throughout  the  flow  field  is  calculated  fret,  the 
Poisson  pressure  equation,  Equation  50,  derived  in  Section  II.  Second 
order  accurate  central  differences  are  used  to  approximate  the  second 
derivatives  of  pressure  in  the  transformed  Laplacian  tern.  The  first 
derivatives  of  pressure  in  the  Laplacian  are  approximated  by  second  order 
accurate  upwind  differences.  The  direction  is  based  on  the  sign  of  the 
coefficients  C  and  t  which  is  the  same  method  used  for  the  momentum  equa¬ 
tions.  In  this  equation,  a  nonlinear  source  term  occurs  which  is  composed 
of  transformation  and  velocity  derivatives.  The  transformation  derivatives 
are  approximated  with  second  order  accurate  central  differences.  Second 
order  accurate  upwind  differences  are  used  for  the  velocity  derivatives. 

The  direction  of  the  differences  at  each  location  is  determined  by  tin- 
direction  of  the  local  \  and  ’  velocity  components  related  to  VC  and  VC 
previously  defined.  The  time  derivative  of  the  velocity  divergence  is 
approximated  with  a  first  order  accurate  backward  difference.  The  require¬ 
ment  that  the  velocity  divergence  at  the  nth  time  ste;  he  zero  is  incor¬ 
porated  into  the  difference  equation  by  setting  Pn  to  zero.  The  numerical 
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non-zero  quantity  D°  ^  is  retained  as  a  correction  factor.  In  this  wav, 
the  static  pressure  is  found  which  tends  to  satisfy  mass  conservation  at 
every  location  for  each  time  step. 

The  following  definitions  are  given  for  the  finite  difference  repre¬ 
sentation  of  the  source  terms  as  described  above: 

I'X  =  (vru  ,  -  y  ait  )/J 

UY  =  (x  jUr  -  xru  {.)/J 

YX  =  (v  v  .  -  v  rv  l/.l 

YY  =  (x  ;vr  -  x,  v  .)  /J 


If  these  definitions,  the  previous  transformation  coefficient  definitions, 
and  index  conventions  are  used,  the  transformed  pressure  equation,  Equa¬ 
tion  V>,  is  rearranged  for  the  field  static  pressure  p? .  given  by 


+  +  2  (YX)  (I’Y)  +  (YY)“  +  AI.FA  (p...  +  n.  P 

L  1+ J  1— 1 


’  ; -*• !  +  Pj-1}  +  BETA(pi+l  j  +  l  "  Pi+1  j-1  +  Pi-3  j-1 

;  - !  is:'  +  rV(4plVl  ■  PIV2)  +  'vv  (4pTYl  -  P.TY2)  '  3(  lT 


*  . )  +  2  (ALFA  +  GAMA ) 


In  this  equation  p..  does  not  appear  implicitly.  The  velocity  component; 


u  j  an!  v  are  only  present  implicitly  in  the  nonlinear  source  tern.s 


The  three  difference  equations.  Equations  92  thru  9a,  form,  a  set  of 
equations  for  the  flow  variables  u,  v,  and  p  at  location  (i,  p  in  tern;;,  of 
the  transformation  derivatives  and  values  of  the  variables  at  neighboring 
points.  The  quasi-linearization  and  differencing  methods  have  decoupled  l h 

ct  to  the  higher  order  terms  for  the  unknowns  at  lor, it  i 
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equations  with 

(  i  ,  j  )  • 


where  s  is  the  iterate  number,  *  designates  the  current  value  obtained 

from  the  difference  equation,  and  e  and  are  acceleration  parameters. 

n  uv  p 

The  acceleration  parameters  are  always  unity  for  the  first  iteration. 

The  computation  of  the  wall  static  pressure  is  carried  out  separately 
as  discussed  in  Section  II.  The  pressure  iteration  technique  of  Cher  in  (hS 
given  by  Equation  40,  avoids  the  difficulty  of  formulating  one  sided  dif¬ 
ferences  at  the  body  for  the  Laplacian  of  pressure.  Also,  at  the  body  sur¬ 
face,  conservation  of  mass  is  directly  imposed.  Second  order  accurate 
upwind  differences  are  used  for  the  transformation  and  velocity  first 
derivatives  in  the  velocity  divergence  of  Equation  51.  The  acceleration 
parameter  1,  which  is  related  to  SOP.  iteration  of  the  Poisson  pressure  oqua 
tion,  is  given  by  Hodge  (53)  as 

I  =  2:..  ,/  At (ALFA  +  GAMA)  ('"■) 

pb 

where  v  is  an  acceleration  parameter  for  the  corresponding  SOL  iteration 
of  the  Poisson  pressure  equation.  For  consistency,  the  iteration  equation 
ir  repeated  as 

(s+i)  =  p.(s)  -  :n  ('-7i 

pb  h 

where  s  represents  tfie  iterate  number  and  the  index  notation  is  kept  as 
described  above. 

The  trailing  edge  pressure  on  an  airfoil  can  lie  difficult  to  obtain 
using  the  iteration  procedure  if  a  sufficiently  fine  grid  is  not  pro: • nt  . 


Hodge  (53)  experienced  such  difficulties  in  his  laminar  flow  work.  The 
trailing  edge  pressure  can  then  be  calculated  using  an  average  01  value- 
obtained  by  linearly  extrapolating  the  pressure;  calculated  at  the  closest 
points  on  the  upper  and  lower  surfaces.  For  a  two  point  linear  extrapola¬ 
tion  in  the  x  direction,  the  pressure  at  the  trailing  edge  becor.es- 


(xl  "  :V  p2  "  (X1  '  V  !'.i 

2(x2  -  x3) 

+  (_^i_  imah-  _PAi '-ax- i_  “  “_xj  mam  - 1 1  1  ’  i : 

"M  x  -  >;  1 

ika:-:-i  imam- : 


(G7a) 


where  the  x's  are  the  x  values;  at  the  selected  surface  points  and  the  sur¬ 
face  j  =  1  value-  is  understood. 

The  iteration  procedure  follows  a  prescribed  order  at  each  point  loca¬ 
tion  (i,  j).  Only  the  body  surface  pressure  is  computed  on  body  points. 

At  field  points,  the  u  velocity  component  is  computed  first,  then  the  v 
velocity  component,  and  finally  the  fi^ld  pressure  p.  The  points  are  order... 
by  varying  the  r  ( j )  index  from  j  =  1  thru  j  -  JMAM-1  or.  successive  ’  I  ir.e- 
in  the  transformed  plane.  The  ;  1  ines  i  =  1  ti.ru  i  =  TMAM-1  ar.  t  raw  r  . 
Values  for  the  variables  at  i  =  1  are  set  equal  to  the  values  at  i  -  I  MAM 
to  insure  continuity  across  the  transformation  branch  cut.  The  i  index 
is  adjusted  across  the  branch  cut  in  a  manner  identical  with  that  described 
for  obtaining  the  transformation  from  equations  83  and  '■>'  . 


C.  Bour’nry  and  Initial  (hpndi  ti  or-. - 

Boundary  and  initial  conditions  are  formulated  for  use  in  the  computa¬ 
tional  plane  in  order  to  completely  define  the  numerical  probler.  Boundar-. 
conditions  on  the  body  surface  ar.  convenient lv  specified  with  tin  bid  - 
fitted  coordinate  svstiT. .  1  he  f  i  r  t  '  contour  is  const  rue  ted  t  .■  *  .  t  he 
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body  surface.  The  analytical  no  slip  boundary  conditions  for  velocity- 
given  by  Equation  52  can  be  directly  applied  to  the  body  surface  in  the 
computational  plane  and  are  given  by 

u  =  0  v  =  0  (98) 

b  b 

where  the  subscript  b  denotes  the  body  surface  y  contour.  The  surface 
pressure  is  calculated  as  discussed  in  Section  IV. B. 

The  r  lines  i  =  1  and  i  =  I MAX  in  the  computational  plane  define  the 
branch  cut  in  the  physical  plane.  No  boundary  conditions  are  required  or 
allowed  on  these  contours.  The  matching  of  solutions  on  these  lines,  which 
was  described  previously,  provides  a  continuous  solution  across  the  cut. 

The  outer  boundary  in  the  physical  plane  frequently  represents  con¬ 
ditions  at  infinity  (freestream  conditions)  given  by  Equations  53  and  54. 

In  the  body-fitted  coordinate  system,  the  outer  boundary  becomes  the  largest 
line  in  the  computational  plane.  The  physical  cuter  boundary  i  -  i  beset! 
to  be  a  circle  with  a  radius  of  10  airfoil  chords.  The  distribution  of 
points  on  the  outer  boundary  is  readily  specified  on  a  circle;  al-o,  th. 
initial  input  solution  for  obtaining  the  t  ran  s  format,  ion  as  discussed  in 
Section  IV. A  is  simple  to  specify.  Chiu  and  Hodge  (82  i  applied  a  nur.i  rive! 
inviscid  anlaysis  to  a  .Joukowski  airfoil  at  angle  of  attack  using  i r t .  - 
stream  boundary  conditions  for  different  outer  boundaries.  They  found 
that  for  a  10°  angle  of  attack  the  lift  coefficient  differed  ’  .  les?  than 
one  percent  and  the  maximum  suction  pressure  coef f icier. t  by  about  two  pi  r- 
cent  from  the  analytical  result.  Na vie r-S t okes  solutions  using  a  circular 
outer  boundary  with  a  smaller  radius  are  compared  wit'  the  J(:  chord  rad  lu¬ 
re  suit  to  determine  effects  of  boundary  placement . 

With  the  outer  boundary  elected  to  approx imat  t  the  far  fit  Id  I n  <  - 
stream,  the  boundary  conditions  for  Velocity  and  pressure  ati  spo  i : ied . 

V? 


As  pointed  out  by  Roache  (83),  caution  must  be  exercised  in  applying  the 
analytical  conditions.  Equations  53  and  54,  which  are  strictly  valid  only 
in  the  limit  of  large  distances  from  the  body.  If  these  conditions  c re¬ 


used  on  a  boundary  at  a  finite  distance  away,  the  result  may  predict  no 
drag  since  a  wake  cannot  exist  at  the  outer  boundary.  This  problem,  is 
avoided  by  applying  upstream  and  downstream  boundary  conditions  on  the 
outer  boundary.  The  upstream  boundary'  along  which  conditions  are  fixed 
is  defined  by  the  semicircle  in  the  half  plane  x  *'  0  where  the  x  axis  lies 
along  the  airfoil  chord  with  the  origin  at  the  midchord.  The  downstream 
boundary  for  which  some  variables  are  permitted  to  be  free  is  the  semicircle 
in  the  half  plane  x  >  0. 

The  upstream  boundary  conditions  for  the  incoming  undisturbed  flow 
become 


u.  =  cos  u 

i  JMAX 


v.  =  sm  'jl 

l  JMA.\ 


»i  JMAX  ‘  °  <"> 

where  i  ranges  over  all  T  contours  which  terminate  on  the  upstream  boundary. 

The  downstream  boundary  conditiors  must  allow  a  wake  downstream  of  a 
body  to  pass  through  the  boundary.  The  no  change  boundary  condition  has 
been  widely  used  (53,83)  and  meets  this  requirement.  One  approach  requires 
no  change  of  the  velocity  components  in  the  f reest ream,  direction.  1‘sing 
the  expressions  for  the  directional  derivative  and  gradient  in  the  trans¬ 
formed  plane,  this  boundary  condition  can  be  written  for  the  u  component 
of  velocity  as 

ur  =  u ,  (xr  sin  «  -  y  cos*)/ (x  .  sin  •  -  v  .  cos  0  (in', 


where  a  is  the  geometric  angle  of  attack.  Mow  consider  a  similar  down¬ 
stream  boundary  condition  where  the  no  change  condition  is  imposed  o:t  the 
velocity  components  along  t  lie  downstream  ’  contours.  In  a  similar  manner, 
this  condition  in  the  transformed  plant  becomes 

UT  =0  V.  -  U  <  1 


This  boundary  condition  is  identical  vith  the  freestream  direction  boundary 
condition  of  Equation  100a  when  the  c  contours  are  in  the  freest rear,  direction. 
In  this  case,  tan  <  =  y  /x  and  Equation  100a  becomes  identical  with  Equation.  I'1 
This  case  does  occur  in  the  downstream  wake  region  for  the  transformation 
obtained  in  Section  IV. A.  Small  changes  in  the  velocity  components  are 
found  in  the  far  wake  region  of  bodies.  The  downstream  far  field  invisciu 
flow  field  also  has  negligible  changes  as  seen  in  the  analysis  of  Appendix 
F.  Thus,  the  no  change  conditions  of  Equation  100b  approximate  the  small  far 
field  downstream  variations  in  velocity  while  allowing  the  momentum  defect 
in  a  wake  to  exit  at  the  outer  boundary. 

The  downstream  outer  boundary  condition  for  pressure  must  he  considered 
along  with  the  velocity  conditions.  Hodge  (53)  investigated  the  use  of  a 
similar  no  change  condition  for  pressure,  namely,  p  =  0  coupled  with  the 
velocity  conditions  of  Equation  100b.  His  laminar  solution  developed  pres¬ 
sure  oscillations.  The  more  restrictive  condition  of  specifying  the  f  ree- 
strear.  pressure  in  the  far  wake  was  successfully  used  by  Hodge  (531.  The 
specification  of  pressure  recovery  with  a  freely  developing  velocity  defect 
in  the  far  wake  gives  a  set  of  numerical  boundary  conditions  which 
can  physciallv  represent  the  far  wake  at  a  finite  far  field  downstream 
boundary.  The  downstream  pressure  boundarv  condition  becomes 

PiJMAv=°  (J01) 

The  suitable  description  of  the  physics  at  the  downstream  computa¬ 
tional  boundary  must  be  written  in  finite  difference  form.  The  no  cii.in,:i 
condition  approximated  with  a  second  order  accurate  central  differer.ee 
applied  at  the  midpoint  in  the  computational  plane  between  the  hour da r. 
point  and  first  interior  point  on  a  downstream  "  contour  (  i  index  •  civ. 
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for  the  downstream  boundary  conditions  of  the  velocity  components.  An. 
alternate  forrnuLation  is  obtained  by  approximating  the  no  change  condition 
with  a  second  order  accurate  upwind  difference  applied  at  the  boundary  point 
Then,  the  downstream  boundary  conditions  for  the  velocity  components  become 
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These  boundary  conditions  also  follow  from  a  second  degree  polynor.inal  fit 
with  a  zero  derivative  imposed  at  the  downstream  boundary. 

The  time  dependent  finite  difference  method  given  previously  require- 
a  set  of  initial  values  for  the  velocities  and  pressures  at  all  the  < i  .  j ) 
locations.  This  requirement  parallels  the  result  discussed  for  the  dif¬ 
ferential  equations  in  Section  lll.b.  Two  methods  are  used  to  provide 
initial  conditions.  In  the  first  method,  an  inviscid  solution  with  cir¬ 
culation  for  the  flow  field  at  the  given  angle  of  attack  is  obtained  using 
SOK  iteration  as  described  in  Appendix  0.  This  solution  then  provide.- 
initial  values  for  the  Xavier-Stokes  solution.  The  second  technique  usis 
the  Navier-Stokes  solution  for  one  angle  of  attack  and  rotates  the  velocity 
by  the  change  in  ancle  of  attack.  The  set  of  initial  values  for  vi.lt  ,  it'  ,->:i 
pressure  is  then  used  to  start  a  Nav  ier-S- okes  solution  at  tin  new  an., -It  of 
a  1 1  a  c  k  . 


D .  Turbul  ence  ‘bnie  l 

The  numerical  procedure  which  implements  the  eddy  viscosity  turbulence 
model  described  in  Section  III  is  discussed.  The  two  lawr  model  near  the 
airfoil  surface  is  presented  followed  by  the  far  wake  mod*  1.  The  proc,  <!ur« 
used  in  the  near  wake  is  then  described.  The  eddv  viscosity  distribut  eci 
throughout  the  flow  field  is  calculated  at  tin  beginning  of  cue!  tint 
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SOK  iteration  is  then  used  to  determine  the  Velocities  and  pressures  at  the 


new  time  level.  The  eddy  viscosity  distribution  is  not  changed  during  the 
iteration  process. 

The  inner  viscosity  computation  from  liquation  6"  requires  values  of 
the  tangential  and  normal  velocity  derivatives  and  the  distance  fro::,  the 
surface.  These  velocity  derivatives  are  given  by  liquations  70  and  71  in 
the  transformed  plane.  They  are  evaluated  using  second  order  accurate  central 
differences  for  both  the  transformation  derivatives  and  the  u  and  v  velocity 
component  derivatives  with  respect  to  •  and  •  .  The  inner  eddy  viscosity  is 
identically  zero  at  the  surface  ('=1)  and  need  not  be  calculated  separate’.;-. 

The  corresponding  velocity  expression  which  is  found  in  the  exponent  of  the 
damping  factor,  Equation  62,  must  be  evaluated  at  the  airfoil  surface  for 
each  r  contour  in  regions  of  turbulence.  The  derivatives  with,  respect  t> 
r  are  approximated  by  second  order  accurate  forward  differences  found  in 
Appendix  B,  and  the  derivatives  with  respect  to  •*  are  evaluated  using  second 
order  accurate  central  differences.  The  distance  from  the  surface  t  •  the 
point  (i,  j)  is  measured  along  the  constant  r  line  represented  by  index 
The  freest  rear,  chord  Reynolds  number  appears  throughout  the  model  a.--  a 
result  of  writing  the  equations  in  none!  irons  icnal  form  and  is  an  input  para- 
met  er . 

The  adverse  pressure  gradient  turbulence  model  modification  in  1  he 
trailing  edge  region  given  by  Equation  to  is  implement.  :  next.  The  eric  1: 
of  the  distance  s,  measured  alone  the  airfoil  surface  iron.  th<  trail  in  •  t  do 
.separation  point,  is  determined  by  examining  the  u  velocity  m::  pox.-nt  at  p  1 :  t  - 
on  tin  •  contour  next  to  tin-  urf.o  »  l  'Mu  sean.li  begin  :  '•  •  1 

cor  resp  xidinc  to  tin-  hS  p,  rcent  chord  1  >>  at  ion  and  adva  cm  toward,  t  !.<  air:  :  1 
lead  in.t  chx  alone  the  second  ’  lit.  h»  carl,  t  err.  i  ::a !  cs  wh<  :  t  h.  . 

value  for  ii .  ,  i  -  do  •  ec  t  ed  v:  i  '  inti  it  <  at  t  a  bed  !'  lev.  he  ri  l.:\.:t  iot.  t  .:  t  i 
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f  is  then  evaluated  for  each  applicable  -  line  where-  trailing  edge  separa¬ 
tion  is  detected.  The  inner  eddy  viscosity  is  multiplied  by  the  rc-laxat  icn 
factor  which  adjusts  the  adverse  pressure  gradient  inner  eddy  viscosit;' 
constant  k,. 

After  the  inner  eddy  viscosity  is  computed,  the  limiting  technique 
evaluates  the  computed  value.  The  inner  eddy  viscosity  is  prevented  from 
decreasing  in  the  outward  normal  and  downstream  tangential  flow  directions 
by  comparing  the  value  with  previously  computed  values  as  given  by  the 
following  sequence: 


ij  -  i  j-1 


ij  -  i-1  j 


The  larger  value  is  selected  during  each  comparison.  The  initial  inner 
eddy  viscosity  limiting  distribution  is  calculated  in  the  attached  boundary 
layer  region  with  a  favorable  pressure  gradient  near  the  airfoil  leading  eth 
T'iie  first  or  outward  direction  limiting  condition  is  imposed  here  as  well. 

The  limit  value  of  the  inner  eddy  viscosity  is  then  compared  with  the 
calculated  outer  layer  eddy  viscosity.  Whenever  the  inner  eddy  viscosity 
first  exceeds  the  outer  value,  the  outer  eddy  viscosity  is:  then,  used  for  tin 
remaining  locations  outward  from  the  surface.  This  switch  insures  the  con¬ 
tinuity  of  the  eddy  viscosity  distribution. 

The  computation  of  the  outer  eddy  viscosity  given  by  Equation,  M  re¬ 
quires  values  for  the  local  boundary  layer  edge  Velocity  and  the  d  i  sp  1  aco::.c: 
thickness  defined  by  Equation  64.  These  quantities  arc  calculated  for  <  a  1 
contour  which  crosses  the  turbulent  boundary  layer.  for  a  :  line,  tin  lei, A 

edge  velocitv  u  is  defined  to  he  the  maximum  tangential  velooitv  r.t  a.- :im  d 
e 


relative  to  the  local  surface  tangent  line.  The  tangential  vc 
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by  Equation  6b  ant!  i~  evaluated  as  described  previously  tor  the  inn.t  r  odd- 
viscositv.  The  boundarv  laver  thickness  '  is  tlic-u  tin  distance  ale-:,  th, 


lir.-  f  run  the  surface  to  the  first  point  win  re  the  tangi  nt  in!  v<  !  o,  it 


less  than  0.9^  u  .  The  displacement  thickness  is  calculated  using  trape- 
o 

zoidal  integration  where  the  limits  of  integrat ion  are  determined  hv  the 
computed  boundary  layer  thickness. 

TVie  resulting  eddv  viscosity  distribution  is  next  modified  to  simulate 
the  transition  from  laminar  to  fully  turbulent  flow.  T he  transition  factor 
T ,  given  by  Equation  67,  is  calculated  as  follows .  The  start  inf  location 
for  transition  on  the  airfoil  surface  must  be  specified.  This  local  ion  is 
specified  by  designating  the  first  grid  point  on  the  surface  downstream  as 
well  as  the  distance  to  that  first  grid  point  .'.s.  The  transition  distance 
s  from  the  starting  location  to  a  given  •'  lino  is  calculated  by  adding  .'.s 
to  the  distance  between  the  successive  surface  grid  points  with  known  loca¬ 
tions.  The  relaxation  factor  is  determined  from  the  transition  factor 
Equation  67  and  the  definition  of  .  If  .  =  3/d  is  specified  when  the 
transition  distance  equals  three  fourths  of  the  total  transition  length,  the 
1/.'.  is  given  by 

1/f  =  2.4456/x  (10-1 

where  Xj  is  the  total  assumed  transition  length.  The  eddy  viscosity,  dis¬ 
tribution  computed  previously  on  a  •"  lint  is  multiplied  by  tin  cor  res;  .cm!  inc 
constant  ”  to  obtain  the  revised  distribution  of  turbulent  eddy  vise.-  it;- 
which  includes  the  transition  region. 

The  final  modification  made  to  the  eddy  viscosity  distribut ion.  acre  ■■■■ 
the  upper  surface  turbulent  region  involves  the  decrease  of  turbult  m  ..  in 
the  direction  from  the  boundary  layer  toward  the  far  field.  7  he  intt  rm ; t  it  n 
factor  giver,  by  Equation  6  3  models  this  behavior .  The  previous  !•  cu  1  .u  !  at  ed: 
boundary  layer  thickness  '  and  the  distance  from  the  airfoil  surto.i  t.  t  h.e 
(i,  j)  location  on  a  line  are  used  to  compute  the  in  t  erm  i  t  teticy  1  o  r  t : . .  t 

lot- at  ion.  1  he  previous  values  of  the  outer  eddv  viscosit  v  are  mult  i  p  I  i  ed:  1 
th*  cort'e  -pond,  ing  vahu  of  t  lie  i  nt  erm  i  1 1  ency  factor  to  obtain  the  J  in.  a  ed  '  . 
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viscosity  distribution  for  the  turbulent  region  above  the  airfoil  sur!a<  <. 

.  tii 

at  trie  n  tine  step. 

Hie  procedure  for  computing  the  eddy  viscosity  distribution  on  a 
line  across  the  turbulent  flow  zone  above  the  airfoil  is  repeated  down¬ 
stream  for  each  line  through  the  95  percent  chord  point.  1  he  contour- 
near  tile  trailing  edge  deviate  from  the  local  normal  direction  in  excess 
of  10'  .  For  locations  in  this  snail  region,  the  eddy  viscosity  is  found 
by  using  t’ie  last  calculated  value  upstream  in  a  direction  tangent  to  the 
air  foil  surface. 

I  he  calculation  of  the  far  wake  eddy  viscosity  given  by  Ernst  for.  7.;  is 
similar  to  the  computation  of  the  Outer  layer  eddy  viscosity.  In  the  wake 
region  the  •  contours  are  approximately  orthogonal  to  the  flow  direction  in. 
the  wake.  The  ’  contours  are  thus  conveniently  defined  paths  across  the  wake 
for  use  in  computing  wake  characteristics. 

i'n  an  •  contour,  the  component  of  velocity  in  the-  frees  treat:  direction 
is  computed  for  the  grid  points  spanning  the  wake .  The  minimum,  of  these 
directional  velocity  components  in  the  wake  is  determined.  The  r.axir 
values  of  the  velocity  components  which,  occur  on.  either  side  of  the  minimur 
value  are  found  and  designated  the  upper  and  lower  edge  velocities .  1 ! *  r. 

the  upper  and  lower  edges  of  the  wake  are  defined  as  the  points  where  t1,., 
local  velocity  component  in  the  freestrear.  direction  f  irst  reaches 
cent  of  the  upper  and  lower  edge  velocities,  respectively .  1  he  distune* •  : r. 

the  point  of  minimum  velocity  to  the  upper  and  lower  wak*  odg* s  art  d*  f in*  d 
as  t  he  upper  and  lower  wake  thicknesses,  respectively.  These  dist.us  * 
determine  the  limits  of  integration  when  calculating  the  upp»  r  and  lower  w.do 
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the  average  of  the  upper  and  lower  displacement  thicknesses.  The  far 
wake  eddy  visco  ;ity  is  calculated  from  Equation  73  using  the  thickness 
and  an  average  of  the  upper  and  lower  edge  velocities. 

The  intermit tency  factor  for  each  grid  pc  ini  on  an  '  contour  is  deter 
mined  using  the  average  value  of  the  upper  and  lower  wake  thicknesses  for 
in  Equation  63.  The  local  distance  to  each  point  is  calculated  from  the 
location  of  minimum  velocity  along  the  '  contour. 

The  near  wake  model  is  numerically  implemented  by  considering  the  i:pp 
turbulent  and  lower  laminar  boundary  layers  separately.  The  eddy  vi sec-sit 
distribution,  which  is  calculated  near  the  trailing  edge  of  the  airfoil  at 
the  95  percent  chord  location,  is  extended  into  the  near  wake  along,  a 
direction  parallel  to  the  airfoil  surface  defined  by  the  93  percent  chard 
and  99  percent  chord  grid  point  locations.  The  chord  lint  is  the  line  of 
extension  for  the  boundary  layers  in  the  near  wake.  The  ’  lines  i  =  1  and 
i  =  IMA.'-'.-l  define  -he  interaction  zone  in  the  near  wake.  The  non- inter¬ 
act  ing  di-  tance  WL  in  Equation  74  is  based  on  the  length  of  the  lover 
laminar  boundary  layer  thickness  computed  at  the  c5  pert  nt  chord  location 
Thi.s  boundary  layer  thickness  is  computed  by  the  same  method  used  previous 
the  upper  surface  boundary  layer.  The  distance  VC. ,  along  the  chord 
. inc  from  the  triiling  edge  to  the  location  where  the  far  wake  eddy  vi.— 
comity  model  begins  is  also  calculated  as  a  multiple  of  the  laminar  trail¬ 
ing  edze  h  jndary  layer  thickness.  The  eddy  viscosity  in  the  interaction 
zone  is  computed  with  the  two  calculated  di. tances  and  the  far  wake  eddy 
viscosity  using  Equation  74. 

E  .  F ore e  Coe f  f  i  <•  i  en  t  s 

The  forces  and  moments  exerted  on  the  airfoil  bode  hv  tin  flow  field 
are  evaluated  at  the  bode  surface  using  Equations  7°,  8m,  and  8.’.  '.hi 


integrals  are  evaluated  using  trapezoidal  integration.  The  product  ter: 
in  the  integrands  of  the  form  (fgl  are  integrated  by  using  products  o: 
averages,  (f.  +  f.  ,1  (g.  -  g.  ,  )/4.  The  r  derivatives  are  evaluated  using 

second  order  accurate  central  differences,  and  the  r  derivatives  are  form 
lated  with  second  order  accurate  upwind  forward  difference?.  The  lift  and 
drag  coefficients  are  calculated  by  adding  the  components  of  and  y  di  recti  or. 
force  coefficients  in  the  normal  freest ream  (90  degrees  counterclockwise)  and 
f reestream.  directions,  respectively,  using  Equations  I1. 13  and  I'. 14. 

The  alternate  technique  based  on  the  control  volume  analysis  in 
Appendix  P  for  computing  the  forces  on  the  airfoil  body  is  numerically 
implemented  in  a  similar  manner.  The  force-  coefficients  in  the  and  y 
directions  for  a  two-dimensional  control  volume  defined  by  an  *  3 im-  arc 

given  by  Equations  83  and  84,  respectively.  The  line  integrals  are  again, 
computed  with  trapezoidal  integration  using  products  of  averages.  The  r 
derivatives  are  evaluated  with  second  order  accurate  central  differences. 

The  ’  derivatives  are  evaluated  with  second  order  accurate  central  dif¬ 
ferences  on  the  interior  '  lines.  Second  order  accurate  upwind  differences 
are  used  to  approximate  the  ’  derivatives  on  the  body  and  outer  boundary 
r  lines.  The  area  integral  is  also  evaluated  using  trapezoidal  integration. 
The  Jacobian  for  an  area  element  bounded  by  the  i  and  i+1  line-'  and  ;-l 
and  j  r  lines  is  calculated  by  taking  the  average  of  the  Jacobian  values 
previously  computed  for  the  viscous  tern,  in  tiie  line  integrals  at  the  cor¬ 
responding  r  segments  on  the  j-1  and  j  ’  lines.  The  area  integrals  arc 
computed  for  both,  n  and  n-1  t imu  intervals  using  the  previously  computed 
flow  field  values.  The  time  deri  vat  iv<_  is  then  evaluated  using  a  *  i  r :  t 
order  accurate  backward  difference. 
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SECTION  V 


I) I  SC C S_SJ  ON  0_4  RESULTS 

Numerical  solutions  for  two-dimensional  incompressible  turbulent 
viscous  flow  over  a  MAC A  0012  airfoil  section  are  obtained  using  the 
implicit  finite  difference  method  previously  described.  The  MACA  0012 
airfoil  section  is  chosen  because  of  its  widespread  use  as  a  test  case 
in  both,  experimental  and  computational  work  (65,66).  The  solutions  arc 
for  angles  of  attack  of  5  ,  7.5",  9.5',  and  11.5"  at  a  chord  Reynolds 
number  of  170,000.  This  Reynolds  number  is  selected  for  comparison  with 
both  previous  NAC'A  data  (5c)  and  for  data  obtained  as  part  of  this  invest i 
gat  ion  at  the  Air  Force  F.  J .  Seiler  Research  Laboratory  (hi,1. 

The  numerically  generated  body-fitted  coordinate  system  which  was 
used  for  each  solution  is  presented.  The  successive-  ver-relaxat ion 
(SOR)  iteration  numerical  parameters  and  convergence  criteria  along  with 
the  steady  state  convergence  criteria  are  then  discussed.  The  velocity 
and  pressure  fields  and  streamline  contours  obtained  from  the  numerical 
solutions  are  examined.  The  calculated  laminar  separation  bubble  charac¬ 
teristics  are  compared  with  empirical  results.  The  position  of  t  ran.-  it  ion 
to  turbulence  and  the  transition  length  relative  to  the  separation  bubble 
location  are  discussed.  The  computed  airfoil  surface  mean  pressure  distri¬ 
butions  are  compared  with  both  an  invisc id  solution  and  experimental  data. 
The  numerically  predicted  lift  and  drag  forces  exerted  by  the  flow  field  o 
the  airfoil  are  presented  along  with  available  experimental  measurement .  . 
The  effect1-  that  the  placement  and  type  of  imposed  far  field  boundarv 
conditions  have  on  the  numerical  solution  are  discussed.  The  sons  it  ivitv 
of  the  numerical  solution  to  grid  size  is  examined. 

h  < 


The  curvilinear  body-fitted  grid  shown  in  Figure  3  is  numerically 
generated  using  boundary  layer  parameters  for  the  attraction  of  '  con¬ 
tours  towards  the  airfoil  surface  and  linear  interpolation  for  wake- 
resolution  given  by  Appendix  C.  The  Blasius  series  flat  plate  chord 
Reynolds  number  of  200,000,  a  nondimensionalized  boundary  layer  edge 
s  velocitv  of  1.2,  and  seven  grid  points  specified  in  the  boundary  layer 
are  used  in  the  ■  contour  attraction  technique  (7).  The  distribution  o! 
airfoil  surface  points  is  accomplished  with  Equation  26  where  A1  =  0.6, 

A2  =  0.4,  and  A3  =  1.0.  Seventy-one  r  or  surface  points  together  with 
forty-four  r  contours  are  used.  The  transformation  difference-  equation-. 
Equations  85  and  86,  were  then  solved  using  SOR  iteration.  The  optimum 
local  acceleration  parameters  varied  between  0.8  and  1.5  throughout  the  f  low 
field  with  essentially  no  change  in  the  local  values  after  25  iteration  . 

After  350  iterations,  the  maximum  errors  in  the  solution  for  both  and 
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y  were  less  than  ,10  ,  and  occurred  near  the  outer  boundary  above  the 

airfoil.  Additional  ~  lines  were  added  in  the  airfoil  waive  region  using 
the  interpolation  method  described  in  Appendix  C.  F.ight  grid  points 
(four  at  a  time)  were  inserted  on  the  rounded  trailing  edge.  Tabic  0.1 
gives  the  •'  and  ’  indexing  for  the  final  79x44  grid  shown  in  Figure  3 
and  the  equivalent  indices  for  the  71x49  original  grid. 

The  airfoil  is  thus  defined  by  a  total  of  79  points.  Eleven  grid 
points  define  the  first  5  chord  nose  region  and  17  points  are  used  to  de¬ 
fine  the  rounded  trailing  edge.  The  maximum  distance  between  sue  ee- i  vi 
surface  grid  points  is  5  of  the  chord.  The  lines  intersect  the.  sur  1  ;n  <. 
contour  within  10"  of  the  local  normal  direction  exc-.pt  in  the  last  5 
chord  region  near  the  trailing  edge. 

The  SOR  iteration  parameters  and  convergence  criteria  which  are  re  a:  i  r*  d 
bv  the  numerical  implicit  procedure  are  similar  for  each  ancle  at  ta.  '■ 
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solution.  The  acceleration  parameters  had  values  of  1.0  for  the  velocity 
components,  0.9  for  the  surface  pessure,  and  1.10  or  1.15  for  field  pressure 
The  larger  value  of  c.  was  used  at  9.5C  and  11.5°  angles  of  attack.  The 
convergence  criteria  for  the  SOR  iteration  required  that  the  maximum  change 
in  the  relative  magnitudes  of  u,  v,  and  p  for  all  locations  be  less  than 
1.1 .  This  criteria  was  relaxed  to  5'"  for  the  9. 5 6  and  11. 5C  cases.  Three 
iterations  per  time  step  were  the  maximum  necessary  during  most  of  each 
computed  case.  For  the  9.5°  case,  the  approach  to  the  steady  state  solu¬ 
tion  was  re-run  with  -u.  ,  =0.95  and  w  =  1.10  for  two  characteristic  times. 

pb  p 

An  additional  iteration  per  time  step  was  required  to  maintain  the  same 
error  tolerance.  Computed  values  of  lift  and  drag  were  virtually  identical. 

Although  implicit  SOR  methods  have  been  shown  by  linear  stability 
analyses  to  be  unconditionally  stable,  convergence  at  a  given  time  re¬ 
quired  a  snail  time  step  At  <  0.001.  The  analysis  in  Appendix  E  based 
on  linear  matrix  theory  indicates  that  the  diagonal  dominance  of  the  set 
of  finite  difference  equations  given  by  Equations  92  and  93  has  a  defin¬ 
ite  dependence  on  1  t  L  t .  Diagonal  dominance  is  a  necessary  and  sufficient 
condition  for  convergence  (81)  of  constant  coefficient  linear  systems 
using  SOR  iteration.  Although  this  system  of  equations  lias  lower  order 
nonlinearities  and  variable  coefficients,  the  order  of  magnitude  study 
of  Appendix  E  indicates  convergence  if  the  most  stringent  requirement 
At  f  0.0003  is  applied.  A  constant  time  step  equal  to  0.0005  was  used 
for  each  solution.  X'o  convergence  problems  were  encountered. 

The  finite  difference  method  described  in  Section  IV  is  a  time  de¬ 
pendent  technique.  The  solution  is  advanced  in  time  by  t  increments. 

The  numerical  solution  was  considered  to  have  reached  a  steady  state  when 
changes  in  the  computed  values  of  the  lift  and  drag  coefficients  were  k  - 
then  1"  over  one  char.ifleri.st  ic  or  nond  irnens  i  nna  1  time  perk'd.  I  it!,  ■-lead.- 
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state  solution  required  three  to  six  characteristic  time-  periods.  and  fro- 
two  to  four  hours  of  CPU  time  on  a  CDC  CYBER  750  computer. 

The  numerical  solution  of  the  mean  flow  field  near  the  airfoil  si  c- 
tion  at  each  angle  of  attack  is  presented  by  using  streamline  contours, 
velocity  field  vectors,  and  pressure  contours.  The  mean  streaml ir.u  con¬ 
tours  are  shown  in  Figures  4  through  7.  A  small  laminar  separation  bubble 
on  the  upper  surface  with  negligible  trailing  edge  separation  is  seen  in 
Figure  4  for  a  =  5C.  For  an  increased  angle  of  attack  equal  to  7.5",  the 
separation  bubble  has  decreased  in  size  and  moved  forward  toward  the  leading 
edge  of  the  airfoil.  Figure  5  also  shows  the  turbulent  trailing  edge  separative 
region  which  has  also  moved  forward.  The  trailing  edge  separation  region 
has  grown  significantly  and  has  progressed  to  the  quarter  chord  position 
when  the  angle  of  attack  has  reached  11. 5C  .  This  behavior  is  typically 
observed  during  trailing  edge  stall  (72) . 

The  velocity  field  vector  plots  presented  in  Figures  8  through  11 
also  illustrate  the  phenomena  discussed  above.  In  addition,  a  sn  1 1  scale 
clockwise  rotational  motion  is  observed  in  the  rear  portion  of  the  trailing 
edge  separation  region.  A  circulatory  flow  pattern  within  the  laminar  separa¬ 
tion  bubble  is  seen  for  >  =  5°  in  Figure  8.  The  boundary  layer  near  the 
leading  edge  for  each  angle  of  attack  contains  at  least  five  grid  points 
along  each  f  contour  while  approximately  ten  points  resolve  the  lover 
surface  laminar  boundary  layer  near  the  trailing  edge.  This  number  if  point: 
should  be  sufficient  to  resolve  the  primary  features.  The  pressun  , o»- 
tours  near  the  airfoil  in  Figures  12  through  15  clearly  show  expand i e  ■ 
regions  of  low  pressure  above  the  airfoil  and  high  pressure  hi  low  the  air¬ 
foil  as  the  angle  of  attack  increases.  The  large  pressuri  variation:  in 
the  upper  surface  suction  peaks  are  also  observed .  The  pressure  i  jo  lb 
indicate  that  only  small  pressure  variations,  occur  in  the  wake  tor  i  .n  ' 
angle  of  attack.  Tlu  snoot!,  and  continuous  numerical  so  1  u  t  i  o'  lot  ;  h* 


pressure  field  is  in  sharp  contrast  to  the  large  oscillatory  results 
obtained  by  Hodge  (7)  for  laminar  flow. 

For  a  moderately  thick  airfoil  with  a  smooth  surface  in  a  flow 
with  a  low  freestream  turbulence  intensity,  the  lan»nar  boundary  layer 
near  the  leading  edge  can  separate  upon  encountering  a  strong  adverse 
'  pressure  gradient.  Subsequently,  shear  layer  instabilities  car:  cause 
transition  to  turbulence.  The  increased  mixing  may  reattach  the  shear 
layer  and  thereby  form  a  separation  bubble.  The  formation  of  separa¬ 
tion  hubbies  and  their  relationship  with  the  stalling  characteristics 
of  airfoils  at  various  Reynolds  numbers  have  been  investigated  by 
Gault  (85),  Gaster  (86),  Hoad,  et  al  (87),  and  Arena  and  Mueller  (?b) . 

The  laminar  separation  bubble  has  been  observed  experimentally  by  Gregory 
and  O'Reilly  (89)  for  the  KACA  0011  airfoil  over  a  large  range  of 
Reynolds  numbers.  These  experimental  and  empirical  results  will  be  used 
to  compare  with  the  numerical  solutions. 

In  this  investigation,  the  beginning  of  turbulent  transition  on 
the  upper  surface  of  the  airfoil  and  the  transition  length  have  been  ba ed 
on  the  closure  of  the  laminar  separating  shear  layer.  The  surface  mean 
pressure  coefficients  for  a  =  5°  and  7.5r  presented  in  Figures  16  and 
17  clearly  show  the  laminar  bubble  constant  pressure  rev  ion,  downs-  trea:. 
of  the  suction  peak,  followed  by  a  rapid  recovery.  Transition  in  tin 
numerical  solution  occurs  near  the  downstream  side  of  the-  bubble  where 
the  steep  pressure  recovery  begins.  This  phenomenon  has  been  r  p  rtid  bv 
Wallis  (90)  and  Arena  and  Mueller  (88).  Die  chord  locations  tor 
start  of  transition  x  and  full  turbulence  x  relat  ive  to  the  .-ep.irat 
and  rent  tacbment  points  of  the-  bubble,  x  and  x,  respet  tivelv,  are  given 
in  Table  1.  For  the  smaller  angles  of  attack,  rea  1 1  achr.en  t  oc,  urri  1 
prior  to  attaining,  fully  turbulent  flow.  This  behavior  has  been  t  >;- 
per i mental 1y  observed  bv  Arena  and  Mueller  (HR i  in  t  hi  i r  low  Revise 'd: 
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TABLE  I 

Computed  Characteristics 

The  Airfoil  Leading  Edge 

a  (Deg) 

xs 

5.0 

1 

O 

IO 

7.5 

-  0.47 

9.5 

- 

11.5 

_ 

of  Turbulence  Near 
(x  =  -0. 5) 


X 

t 

XR 

XT 

0.24 

-  0.15 

-  O.no 

0.42 

-  0.38 

-  0. 35 

0.48 

- 

-  0.40 

0.49 

_ 

-  0.48 

TABLE  II  Computed  Laminar  Separation  Bubble  Characteristics 


■i  (Deg) 

M 

S 

H12S 

Ref  *s 

Re  e 

'  s 

"R 

5 

.124 

3.1 

6S0 

o  o  q 

.002  5 

.27 

-  .  002 

7.5 

.091 

3.2 

540 

2  68 

.0013 

.09 

-  .002 

number  flow  study. 

The  computed  values  of  several  parameters  at  separation  for  both 
a  =  5°  and  a  =  7.5°  cases  agree  with  empirical  results  and  are  presented 
in  Table  11.  The  pressure  gradient  parameter  >1  =  -  Re  du  /c!s 

l' 

in  nond intensions 1  variables,  where  is  the  momentum  thickness  and  s 

is  arclength,  compares  favorably  with  the  laminar  separation  criteria 

of  Curie  and  Skan  (91)  given  by  M  >_  0.09.  The  shape  parameter  H.,,  = 

£*/.  takes  on  values  of  3.1  and  3.2  at  laminar  separation  compared  with 

an  average  value  of  3.5  reported  by  Curie  and  Skan.  The  Reynolds  number 

at  separation  based  on  the  displacement  thickness  (Re;.j.g)  lias  a  value 

greater  than  500  and  predicts  a  short  bubble  using  the  Owen-Klanfer  (°2) 

criteria  confirmed  by  Gault  (85)  and  Crabtree  (93).  The  bubble  length  1 

is  of  the  order  10“  £*  and  decreases  rapidly  with  increased  angle  of  attack 

as  seen  in  Table  11.  This  relationship  also  indicates  a  short  bubble  (90, 

U 

93)  whereas  long  bubbles  have  lengths  of  the  order  10  ?*  (85,8b).  The 

assumption  that  laminar  separation  precedes  turbulent  transition  is  verified 
with  the  computed  Reynolds  number  at  separation  based  on  momentum  thickness 
Re  .  Crabtree's  (99)  criteria  indicates  that  transition  has  occurred  if 
Re-  j.  >  780  when  M  first  reaches  0.09.  Thus,  the  solutions  correctly  predict 
that  laminar  separation  occurs  before  transition  to  turbulence.  The  turbu¬ 
lent  reattachment  criterion  for  the  bubble  given  by  Roberts  (‘'5'  is  = 

(  I w  )  du  /tls  >  -  0.0059.  This  criterion  is  also  satisfied  as  show:',  in 
e  e  - 

Table  II. 

The  trends  discussed  above  for  M.. ,  and  both  boundary  lav»  r  h.  vnold- 

S  K 

numbers  continue  at  the  higher  angles  of  attack.  However,  at  a::,;l,s  si 
attack  near  stall  the  bubble  length  decrea?  os  to  about  1  chord  <  '<')  >  which 
is  the  ordi  i  of  the  ’  contour  spacing  nea:  the  leading  <  dee.  Th«  m::s  r  i  i  5 
solution:  at  ■  =  9.5  and  <  ~  11.5  consequently  do  not  resolve  t  h>  bid’  1.  ,  .. 


seen  in  Figures  13  and  19.  The  solutions  become  very  sensitive  to  small 
changes  in  the  turbulent  transition  length  and  starting  location  .  hr 
example,  at  a  =  9.5°  with  an  initial  transition  location  of  >:  =  -0.48, 

a  change  in  transition  length  from  0.022  to  0.025  produced  shedding 
of  hubbies,  large  pressure  oscillations,  and  a  significant  30'  increase 
in  average  drag  compared  with  the  steady  state  solution.  The  experimental 
pressure  data  were  measured  with  diaphragm  transducers  capable  of  measuring 
frequencies  well  beyond  5000  Hz.  No  dominant  frequency  was  observed. 

Hoad,  et  al  (87)  recently  observed  a  small  scale  unsteadiness  in  a  short 
region  near  the  surface  of  a  NACA  0012  airfoil  at  angle  of  attack.  The 
laser  velocimeter  mean  velocity  histograms  for  a  certain  region  near  the 
leading  edge  exhibited  dominant  and  small  secondary  mean  velocities.  Head, 
et  al  suggest  that  this  behavior  may  indicate  an  unsteadiness  within  the 
laminar  separation  bubble.  The  remaining  flow  field,  however,  was  steady. 

The  numerical  solutions  presented  for  each  angle  of  attack  are  the 
result  of  a  parametric  study  where  the  upper  surface  transition  length  and 
location  were  varied  in  order  to  close  the  bubble  and  maintain  steady  flow. 
In  each  case,  a  further  increase  in  the  transition  length  or  movement  down¬ 
stream  of  the  start  of  transition  resulted  in  a  non-physical  oscillator; 
motion  propagated  downstream  with  an  accompanied  large  increase  in  dm.-. 

Tiie  good  agreement  between  the  calculated  and  empirical  separation  hid!  1 , 
characteristics  indicates  that  the  origin  of  turbulence  or.  the  upper  sur¬ 
face  is  satisfactorily  modelled. 

Computed  boundary  layer  veloi  itv  profiles  at  four  stations  on.  t  ■ 
upper  surface  are  compared  with  two  sets  of  hot  wire  anemone t rv  data  (.■■'•'  is 
Figures  20  through  23.  The  profiles  were  measured  and  computid  at  tin  an  < 
locations  on  constant  lines  located  approx imat  el v  at  It  ,  Vs  ,  f  '•  ,  and 
°2  chord.  The  profiles  art  nond  i  mens  i  on.i  !  i  fed  1-v  the  locally  co-  put,  1 
houndurv  Inver  thickms  which  is  defined  a  .  th.  noma!  .li'tas.i  :  rot-  tin 
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surface  to  where  the  tangential  velocity  f  ir.-t  attains  u'-j  of  tie. 
edge  velocitv.  Significant  RMS  fluctuations  for  each  angle  <  :  at  t  ,if .  w.  t. 
detected  experimentally  at  the  inner  data  point  location.-  5  or  •  tat  i.-:>  tw. 
through  four  which  indicate  the  presence  of  a  turbulent  boundary  l.r-i 
No  fluctuations  were  observed  in  the  lower  suri  ao«.  hcundar;  l.v. .  r  .w. 
the  90'  chord  location.  This  result  provides  experimental  t  vid.  i  :  •• 
the  assumption  that  the  lower  surface  boundary  layer  remains  laminar  it: 
the  presence  of  the  favorable  pressure  gradient  at  the  -Ids  num!  t  r 

under  investigation.  The  calculated  boundary  layer  becomes  third,  r  down¬ 
stream  compared  with  the  experimental  profiles.  ini,-,  result  could  In. 
caused  by  grid  boundary  layer  resolution  errors  propagated  down.-:  rear  .  r 
by  deficiencies  in  the  tubulence  model.  The  hot  wire  data  i.avi  an  et.- 
mated  experimental  error  of  5' . 

The  boundary  layers  on  the  upper  and  lower  .surfaces,  merge  to  for::,  t  hi 
wake  downstream  of  the  airfoil.  Wake  meat:  velo. it y  profiles  at  chord  lo¬ 
cations  of  0.58,  0.79,  and  1.54  for  each  ancle  of  attack  art  shown  in. 
Figures  24  through  27.  The  trailing  edge  is  located  at  “  !  .5.  ihe  ;>v.  - 
files  are  measured  along  constant  r  lines  with  tin  origin  on  the  extend.,  d 
chordline.  The  thicker  numerical  upper  surface  boundary  layer  propagates 
downstream.  The  computed  large  variations  of  velocity  near  the  low. r  u;a 
of  the-  wake,  which  come  from  the  lower  surface  boundary  layer,  are  in  go.  d 
agreement  with  the  hot  wire  anenometry  measurements.  The  measurement--  wv  u 
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1  A:;:.!  II!  Computed  Far  Wake  Charaet  er  ist  ii  •• 
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■}.,  he -'acids  stress  u  v  ,  nontira-ns  iun.i  1  i  neU  :  ,  was,  v  .  :  t  iJ 

and  contour  plot-  arc  shown  in  Figures  28  through  31.  1  hi-  geor.etr  i, 

p  a  1 1  e  r  r.  qualitative!.'  si~il.ar  to  t  ne  oxpcr  inie.nta  1  data  obtained  :  y 

Coles  and  W'adcoek  (9(1  lor  a  NACA  4412  airfoil  which  arc-  shown  in 
i  i pure  72.  Some  experimental  measurements  were  made  during  the  hot 
wire  anemor.o t ry  study  (84)  and  are’  presented  in  Appendix  "his  da:  a 

also  illustrates  the  rapid  variation  of  negative  and  positive  strossi - 
in  the  upper  and  lower  portions  of  the  near  wake,  respectively .  -v 
tudes  of  the  order  0.01  have  been  observed  (84,  961  in  the  near  wake  and 
were  also  obtained  numerically  where  the  contour  values  vary  from  — !  . r !  1  U 
+  ’.0l  in  Figures  '28  through  71. 

The  computed  surface  mean  pressure  distributions  for  all  four  angU  ■ 
of  attack,  presented  in  Figures  ] 6  through  1°,  favorably  compare  wit:, 
the  experimental  data  (84)  which  have  an  estimated  error  of  7  to  7  .  A", 

invise  id  solution  wit!’  imposed  Kutta  condition  was  computed  using  t:.e 
approach  outlined  in  ‘.ppendix  0  and  is  also  shown.  The  cor..;  ..led  hava  r- 
S tokos  laminar  suction  peak  is  well  defined  and  increases  rapid iv  wit:, 
ancle  of  at'  ack .  The  lower  surface  pressure  peak  hot!,  moves  dovro  treat: 
and  become--  broader  as  the  angle  of  attack  increases.  1’  is  h  :,:vc  r  i . 
also  seen  in  the  experimental  data.  A  separation  Hdd  K  is  o'  m  ro  in. 
the  expe  r  i  nu  n  *  a  1  data  where  the  curvature  rapidly  changes  .  • :  •  hut:  ■ . 

occurs  further  d  own  st  rear  when  compared  with  the  N  rcicr-ht  >  -ki  1  ;  ut  at  ;  .»• 

Tli  i  --  d  i  f ‘'erence  is  probably  caused  by  free  tn.r  t  urhul  era.  c  w!.e!.  1  ■■’:  d,. 
s.  "..J  r.j :  lin;  and  t  i . i  -  si;'i  .i  pier.'.  !  C  I"  a  t  i  o:'  ' f  tie-  bn1'!  .  >  .  •  '■  i  '  t  .  •  ’ : :  .i  i 


T 


data  were  measured  in  a  wind  tunnel  with  a  fret-stream  turbulence  intern.- il 
o:  approximately  0.5  compared  with  the  unperturbed  freest  ream  of  the  run  - 
erica  1  solution.  Some  disagreement  occurs  in  the  trailing  edge  region 
where  the  computed  result  has  a  more  pronounced  flat  trailing  edge  stall 
character  ist  ic . 

The  computed  inviscid  surface  pressure  distribution  consistent!'' 
predicts  both  a  larger  suction  peak  on  the  upper  surface  and  a  larger  pres¬ 
sure  peak  on  the  lower  surface.  The  discrepancies  increase  wit?',  angle  of 
attack.  The  inviscid  solutions  predict  complete  recovery  at  the  trailing 
edge  with  a  stagnation  point  and  C  =  1.0.  The  large  pressure-  gradier.t.- 
in  the  inviscid  and  viscous  solutions  near  the  airfoil  nose  are  resolved 
which  indicates  an  adequate  grid  point  distribution  in  the  region. 

An  examination  of  the  experimental  and  inviscid  surface  pressure 
distributions  for  zero  angle  of  attack  given  in  Figure  33  provides 
additional  information.  The  inviscid  pressure  distributions,  obtained 
with  the  S:'*K  finite  difference  technique  in  Appendix  C,  for  the  upper 
and  lower  surfaces  are  indistinguishable.  This  physically  correct 
symmetric  result  came-  direct.;.-  iron  the  numerical  solution,  and  was  not 
specified  as  a  boundary  condition.  Furthermore,  the  numerical  inviscid 
result  is  virtually  identical  with  an  inviscid  solution  ( •  71  using  t  he 
method  of  Theodorsen  (12).  The  experimental  data  for  both  surfaces, 
agree  well  with  the  inviscid  results  which  indicates  the  small  offer  t 
that  viscosity  has  on  the  pressure  d  istribnt  ion  for  a  stream!  in.ed  ■--  - 

metric  airfoil  at  zero  angle  of  attack.  This  agreement  a  1 provide.-, 
evidence  that  the  experimental  airfoil  sect  ion  does  ap;  -oxira: t  t h.  FA 
On  1  2  ror.i  igntatiun.  Viscous  i  fleet--  on  the  pressure  diet  r i  1  v ?  Ion  do  -  .  .  nr 
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reveals  a  slight  suction  peak  on  the  lover  surface .  1  he  suction  peak 

suggests  that  the  experimental  zero  ang  1  e  of  attack  n.ay  actually  be 
slightly  negative.  A  smal  deviation  in  angle  is  plausil le  since  the 
experimental  zero  was  determined  by  manually  adjusting  the  airfoil 
attitude  until  the  pressure  distributions  on  loth  surfaces  appear*.  - 
idea  l  i  ca  1  . 

The  lift  coefficients  obtained  from  the  present  turbulent  h'avii  r- 
Stokes  solution,  inviscid  results,  and  two  sets  of  experimental  data  ar* 
compared  in  Figure  is.  The  experimental  results  shown  have  been  icrre.h 
for  angle  of  attack.  The  effective  experimental  zero  ancle  of  attack  for 
a  symmetric  airfoil  section  is  defined  to  he  the  angle  where  the  lift 
coefficient  is  zero.  For  the  Seiler  Laboratory  results  th*  effec¬ 

tive  zero  was  found  to  he  0.5s.  This  snail  correction  is  in  accord  with 
the  analysis  of  the  nominal  zero  angle  of  attack  pressure  distribution. 


results  already  presented. 


xperimer.tal  data  have  therefore  V 
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results  ( S b  were  sir.ilarlv  published  in  corrected  fort 
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two-dimensional  wind  tunnel  corrections  including  solid  Mocking,  wake 
Mocking,  and  streamline  curvature  effects  We-r*  invest  igated.  Fn; irical 
low-speed  wind  tunnel  correction  factors  (97‘>  using  the  M  Her  1  h:  rat  v 
d'xV  tunnel  geometry  were  applied  to  the  experimental  data  at  II  H 
account  for  blockage  and  wall  int  erf  <  ren.ee.  The  cor  re  .'ted  value 
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The  NACA  lift  curve  then  flattens  more  rapidly  near  stall.  These  minor 
differences  may  be  explained  by  a  larger  freestream  turbulence  level  in 
the  NACA  results  which  tends  to  delay  the  turbulent  trailing  edge  separa¬ 
tion.  Also,  the  major  variation  of  the  lift  curve  at  small  angles  of 
attack  occurs  for  Reynolds  numbers  less  than  10^  (56,  89).  Freestream 
turbulence  can  cause  more  rapid  transition  and  result  in  an  apparent 
freestream  Reynolds  number  greater  than  the  nominal  turbulence-free  value. 

The  computed  turbulent  Navier-Stokes  lift  coefficients  for  the 
four  angle  of  attack  cases  are  in  excellent  agreement  with  the  experi¬ 
mental  data  as  shown  in  Figure  34.  The  force  components  were  computed 
using  trapezoidal  integration  to  sum  the  total  surface  stresses  obtained 
from  the  flow  field  solution.  The  numerical  results  exhibit  a  gradual 
decrease  in  slope  with  increased  angle  of  attack  similar  to  the  Seiler 
Laboratory  data.  This  behavior  is  consistent  with  the  data  because  the 
numerical  solution  has  a  quiescent  incoming  freestream  with  a  correspond¬ 
ing  smaller  apparent  Reynolds  number.  The  significant  effect  that  the 
viscous  separation  in  the  trailing  edge  region  has  on  the  lift  coefficient 
for  angles  of  attack  greater  than  8°  is  observed.  The  computed  values 
are  within  5%  of  the  experimental  results.  The  lift  coefficient  was  also 
computed  for  each  angle  of  attack  using  the  contour  momentum  integral 
method  described  in  Appendix  D.  The  calculated  values  differ  by  less 
than  1%  for  n  contour  paths  within  one-half  chord  of  the  aifoil  surface. 
This  result  demonstrates  the  near  flow  field  consistency  in  resolving  the 
lift  force.  The  time  dependent  terms  in  Equations  83  and  84  always  con¬ 
tributed  less  than  0.5%  of  the  total  computed  lift  w’hich  indicates  again 
a  converged  steady  solution. 

Two  inviscid  flow  predictions  for  the  lift  coefficient  are  shown  in 
Figure  34  for  comparison.  The  numerical  inviscid  finite  difference  results 
were  computed  with  the  second  order  accurate  trapezoidal  integration 
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technique  used  in  the  viscous  calculations.  The  linear  airfoil  theory 
2tt  slope  is  included.  The  inviscid  flow  results  over-estimate  lift  by 
25%  at  an  angle  of  attack  of  5°.  This  discrepancy  increases  with  angle 
of  attack  because  the  viscous  effects  which  induce  stall  are  not  pre¬ 
sent  in  an  inviscid  calculation. 

The  computed  drag  coefficients  for  the  four  cases  are  compared  with 
available  experimental  data  (56)  in  Figure  35.  The  force  components  were 
computed  using  trapezoidal  integration  to  sum  the  total  surface  stresses 
obtained  from  the  flow  field  solution.  The  drag  component  is  in  the 
direction  of  the  incoming  freestream  flow.  The  rapid  increase  in  the 
drag  coefficient  which  accompanies  the  smaller  increase  in  lift  at  higher 
angles  of  attack  near  stall  is  observed.  Jacobs  and  Sherman  (56)  acknow¬ 
ledged  that  the  experimental  drag  data  in  this  Reynolds  number  region 
have  error  greater  than  the  estimated  +  0.001  for  the  larger  Reynolds 
number  results.  The  agreement  between  the  present  numerical  solution  and 
the  experimental  data  is  within  ten  drag  counts  in  the  region  of  the 
maximum  lift  to  drag  ratio. 

The  presented  numerical  solutions  for  two-dimensional  incompressible 
turbulent  viscous  flow  over  airfoils  were  obtained  with  several  parametric 
studies.  The  effects  associated  with  varying  the  turbulence  transition 
values  and  the  SOR  iteration  parameters  including  time  step  size  have 
already  been  discussed.  The  influence  of  other  turbulence  model  para¬ 
meters,  far  field  boundary  conditions,  and  grid  fineness  on  the  numerical 
solutions  have  also  been  investigated. 

The  turbulence  model  contains  several  parameters  which  can  vary.  The 
effects  that  changes  in  these  parameters  have  on  the  flow  field  near  the  air¬ 
foil  surface  and  on  the  values  of  lift  and  drag  were  examined.  The  value 
of  the  nominal  inner  eddy  viscosity  parameter  k^  was  increased  by  20°'  at 
a  5°  angle  of  attack.  The  lift  coefficient  subsequently  increased  by  2°/ 
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while  the  drag  coefficient  increased  by  57.  or  six  drag  counts.  These 
integrated  effects  are  caused  by  the  observed  small  increase  of  velocity 
in  the  inner  boundary  layer  on  the  upper  surface.  The  outer  eddy  vis¬ 
cosity  parameter  was  next  increased  by  50%.  The  result  was  an  increase 
in  turbulent  mixing  with  a  corresponding  thicker  boundary  layer.  The 
velocities  near  the  edge  of  the  upper  surface  boundary  layer  decreased 
which  caused  an  8%  loss  of  lift  with  no  apparent  effect  on  the  calculated 
drag. 

Values  for  the  constants  in  the  inner  eddy  viscosity  parameter  k^ 
relaxation  factor  given  by  Equation  65  were  obtained  from  a  parametric 
study  at  an  angle  of  attack  of  11.5°.  This  angle  of  attack  was  chosen 
because  the  effects  of  the  relaxation  are  significant  only  for  angles  of 
attack  near  stall.  A  solution  was  initially  computed  without  the  relaxa¬ 
tion  factor  f.  The  laminar  separation  bubble  remained  closed  with 
trailing  edge  separation  beginning  at  mid-chord.  The  calculated  values 
for  the  lift  and  drag  coefficients  were  0.98  and  0.064,  respectively.  An 
examination  of  the  lift  curve  in  Figure  34  reveals  that  this  solution  ex¬ 
hibits  the  character  ox  leading  edge  stall  where  an  approximate  linear 
behavior  is  sustained  until  separation  abruptly  occurs.  Leading  edge  stall 
does  occur  for  the  NACA  0012  airfoil  for  Reynolds  numbers  greater  than 
about  500,000  (56,  89).  Thus  the  mechanism  for  leading  edge  stall  seems 
present  within  the  modified  turbulence  model.  A  numerical  solution  was 
next  obtained  using  the  relaxation  factor  with  a  relaxation  distance  = 
0.25  and  a  delay  distance  s^  =  0.  in  nondimens ional  distances.  The  laminar 
bubble  remain  closed,  but  the  trailing  edge  separation  region  was  signifi¬ 
cantly  larger  and  moved  to  within  the  20%  chord  location  (x  =  -  0.3).  The 
lift  coefficient  decreased  substantially  to  a  value  of  0.8  while  the  drag 
coefficient  increased  moderately  to  a  value  of  0.074.  The  final  solution 
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which  has  been  previously  discussed  was  calculated  using  distances  = 

0.5  and  s^  =  0.1.  The  computed  lift  and  drag  coefficients  became  0.84 
and  0.068,  respectively,  and  are  in  agreement  with  the  experimental  data. 
The  sensitivity  of  the  numerical  solution  to  these  parameters  has  thus 
been  obtained.  The  last  set  of  distance  parameters  was  used  in  calcu¬ 
lating  the  solutions  for  the  other  angles  of  attack.  The  onset  of 
trailing  edge  separation  for  the  5°  and  7.5°  solutions,  seen  in  Figures 
4  and  5,  indicates  that  the  relaxation  factor  was  not  activated.  Thus 
the  relaxation  model  for  parameter  affects  the  flow  field  only  near 
stall  where  the  pressure  gradient  approaches  zero  over  a  large  portion 
of  the  upper  surface. 

The  relaxation  distance  WL^  in  the  near  wake  turbulence  model  given 
by  Equation  74  was  also  investigated.  Numerical  solutions  were  obtained 
at  5°  angle  of  attack  for  relaxation  distances  equal  to  5  and  100  times 
the  lower  surface  boundary  layer  thickness  near  the  trailing  edge.  No 
changes  in  either  the  surface  pressure  distribution  or  the  integrated 
force  coefficients  were  observed.  In  two  recent  near  wake  calculations 
with  similar  relaxation  models,  Waskiewicz  (98)  used  a  value  of  30  trail¬ 
ing  edge  boundary  layer  thicknesses  and  Hasen  (99)  used  10  boundary  layer 
thicknesses.  The  relaxation  distance  of  5  boundary  layer  thicknesses  was 
retained  since  the  agreement  of  the  Reynolds  stress  field  discussed  pre¬ 
viously  was  improved  for  this  choice. 

The  effect  on  the  solution  of  the  placement  and  type  of  far  field 
boundary  conditions  should  be  considered  in  any  computational  work.  In 
this  investigation,  two  types  of  boundary  conditions  at  different  loca¬ 
tions  were  examined.  The  four  angle  of  attack  solutions  that  have  been 
discussed  were  computed  using  the  freestream  far  field  boundary  conditions 
described  in  Sections  III.B  and  IV. C  at  a  circular  outer  boundary  of  radius 
10  chords.  The  effect  of  the  far  field  boundary  placement  on  an  inviscid 


solution  for  Joukowski  airfoils  was  previously  accomplished  (82) . 

This  analysis  indicated  that  a  radius  of  10  chords  was  sufficient  as 
discussed  in  Section  IV. C.  The  effect  of  the  boundary  placement  on 
the  present  viscous  solution  was  examined  by  using  an  alternate  outer 
boundary.  The  J  =  40  0  contour  with  a  semi-major  x  axis  of  4.77  and 
semi-minor  y  axis  of  4.64  was  chosen  as  the  new  outer  boundary  because 
it  approximates  a  circle  with  half  the  previous  radius.  The  numerical 
solution  was  then  obtained  for  a  =  5°  and  compared  with  the  solution 
using  the  10  chord  radius  outer  boundary.  The  average  values  for  the 
lift  and  drag  coefficients  were  identical.  The  solution  with  the  closer 
far  field  boundary  had  a  small  oscillatory  behavior  with  variations  in 
the  lift  and  drag  coefficients  of  1%  and  2%,  respectively.  Further  ex¬ 
amination  of  the  near  surface  flow  field  revealed  that  the  laminar 
separation  bubble  had  a  small  scale  unsteadiness  which  locally  affected 
the  pressure  field. 

The  far  field  potential  flow  boundary  condition  model  developed  in 
Appendix  F  was  next  applied  with  the  outer  boundary  locations  of  10  chords 
and  about  5  chords  in  turn.  This  model  approximates  the  small  pertur¬ 
bations  from  the  freestream  conditions  for  the  velocity  and  pressure 
fields  at  a  large  but  finite  distance  caused  by  the  presence  of  a  body  in 
the  flow  field.  In  this  approach,  the  upstream  boundary  conditions  for 
both  velocity  and  pressure  became  the  calculated  values  from  the  inviscid 
potential  flow  model.  The  downstream  boundary  condition  for  pressure 
also  became  the  corrected  value  rather  than  the  freesteam  value.  The  no 
change  downstream  boundary  condition  for  the  velocity  components  was  re¬ 
tained.  The  solution  was  calculated  with  the  revised  far  field  boundary 
conditions  at  the  10  chord  radius  circular  outer  boundary.  The  boundary 
conditions  were  initially  updated  every  200  time  steps  (0.2T)  by  using 
the  latest  calculated  value  of  the  circulation.  After  three  characteristic 
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time  periods  T,  the  mean  values  of  the  lift  and  drag  coefficients 
were  virtually  identical  with  the  previous  freestream  boundary 
condition  results.  However,  while  the  lift  coefficient  had  variations 
of  less  than  1%  from  the  mean  value,  the  calculated  drag  coefficient  was 
periodic  with  a  period  of  0.4T  and  variations  of  +  10".  Since  the  period 
of  the  oscillations  was  twice  the  period  of  the  boundary  condition  update 
and  the  circulation  was  essentially  a  constant  at  this  time,  the  outer 
boundary  conditions  (except  the  no  change  downstream  condition)  were  kept 
at  the  current  values  and  the  solution  was  advanced  1.5  characteristic 
time  periods.  At  this  time,  the  lift  coefficient  had  a  steady  value  of 
0.425  and  the  drag  coefficient  became  0.011  +  2%  (two  drag  counts). 

These  computed  results  are  within  12  of  the  corresponding  coefficients 
calculated  from  the  solution  using  the  simple  freestream  conditions. 

Thus  the  more  accurate  boundary  condition  which  includes  the  effect  of 
a  finite  distance  from  the  body  induces  very  small  changes  in  the  computed 
flow  field  near  the  airfoil. 

The  modified  far  field  potential  flow  boundary  conditions  were  also 
used  at  the  near-circular  five  chord  radius  outer  boundary  defined  above 
to  obtain  a  solution  again  for  Q  =  5°.  After  two  characteristic  times, 
the  outer  modified  freestream  boundary  conditions  were  held  constant  since 
the  lift  coefficient  had  small  variations  less  than  1%.  The  solution  ex¬ 
hibited  a  small  oscillatory  behavior  similar  to  the  freestream  boundary 
condition  result  even  after  four  characteristic  time  periods.  The  lift 
coefficient  had  a  mean  value  within  1%  of  the  freestream  boundary  condition 
value  with  variations  below  1%.  The  drag  coefficient  had  the  same  average 
value  as  the  previous  result,  but  the  oscillations  were  larger  with  a  de¬ 
viation  from  the  mean  of  approximately  +  10".  (10  drag  counts).  The  un¬ 
steadiness  was  again  observed  to  emanate  from  the  laminar  separation  bubble. 


Small  pressure  oscillations  propagated  downstream  on  the  upper  surface 
of  the  airfoil  and  were  eventually  damped  out  near  the  trailing  edge.  In 
eac i  case  no  changes  were  made  in  any  of  the  turbulence  or  other  parameters 
during  the  computations.  Small  changes  in  the  turbulence  transition  near 
the  bubble  would  probably  prevent  the  unsteadiness. 

The  far  field  boundary  condition  study  has  indicated  that  the  use 
of  the  freestream  boundary  conditions  at  the  large  but  finite  distance 
from  the  airfoil  surface  is  sufficient  for  the  computation  of  the  near 
surface  flow  field  and  force  coefficients.  Variations  in  the  outer  bound¬ 
ary  placement  and  the  use  of  a  more  accurate  far  field  boundary  condition 
had  negligible  effects  on  the  present  numerical  solution. 

The  numerical  implementation  of  the  downstream  no  change  boundary 
condition  was  also  investigated.  Second  order  accurate  central  spatial 
and  upwind  differences  were  used  in  separate  solutions  for  an  angle  of 
attack  of  7.5°.  No  difference  (+0.01/0  was  observed  in  the  surface  pres¬ 
sure  distributions  or  computed  force  coefficients  between  the  two  bound¬ 
ary  conditions.  The  computed  velocities  at  locations  across  the  wake  on 
the  outer  boundary  for  the  two  cases  differed  by  less  than  17.  The  upwind 
difference  formulation  was  chosen  because  of  the  associated  convective 
properties . 

The  sensitivity  of  the  numerical  solution  to  the  fineness  of  the 
grid  was  examined.  The  7.5°  angle  of  attack  case  was  chosen  because 
both  the  laminar  separation  bubble  and  a  region  of  separated  flow  near 
the  trailing  edge  are  present.  A  coarse  grid  was  obtained  from  the  79x44 
grid  by  deleting  the  odd  numbered  n  contours  except  for  the  body  contour. 
The  coarse  79x23  grid  was  a  subset  of  the  full  coordinate  system  and  was 
used  to  study  the  effect  of  boundary  layer  resolution  on  the  numerical 
solution.  The  turbulence  parameters  3nd  SOR  acceleration  parameters  were 
unchanged  from  the  previous  solution.  The  coarse  grid  computation  was 
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carried  out  over  four  chacteristic  time  periods.  The  calculated  values 
for  the  lift  and  drag  coefficients  approached  a  "steady"  periodic  state 
after  two  characterise  times  with  a  period  of  0.6T.  The  mean  value  of  the 
lift  coefficient,  averaged  over  the  last  two  characteristic  times,  became 
0.606  +  27  compared  with  the  full  grid  solution  value  of  0.606.  The 
computed  mean  drag  coefficient  was  0.0260  +  15%  compared  with  the  79x44 
grid  solution  result  of  0.0240.  Values  for  the  wake  velocities  in  the 
freestream  direction  along  the  outer  boundary  for  each  solution  compared 
within  1%.  The  unsteadiness  was  found  to  originate  from  the  oscillating 
laminar  separating  shear  layer  which  forms  the  front  portion  of  the  bubble. 
The  edge  of  the  boundary  layer  varied  between  the  J=3  and  J=4  0  contours, 
as  defined  by  the  79x23  grid  system,  at  the  1%  chord  location.  This  motion 
also  caused  the  bubble  to  vary  in  length  and  produced  pressure  oscillations 
along  the  upper  surface.  This  unsteadiness  is  probably  attributed  to  the 
reduced  resolution  of  the  bubble  fn  the  0  direction.  Similar  sensitivity 
has  already  been  discussed  concerning  the  resolution  of  the  extremely 
short  bubble  in  the  ?  direction  at  larger  angles  of  attack. 

The  computed  flow  field  near  the  airfoil  and  the  calculation  of  the 
force  coefficients  using  the  coarse  grid  are  in  agreement  with  the  results 
for  the  complete  79x44  grid.  A  quadratic  Richardson's  extrapolation  on  the 
computed  drag  coefficients  indicates  an  error  of  +0.0007.  This  agreement 
indicates  that  the  solution  obtained  with  the  full  79x44  coordinate  system 
is  adequate  to  obtain  +10  drag  counts  and  therefore  sufficiently  approximates 
the  limiting  numerical  solution  obtained  from  an  extemely  fine  grid.  The 
coarse  grid  computation  further  indicates  that  the  numerical  multi-grid 
approach  (100)  may  be  implemented  for  complex  flows  using  general  grid  trans¬ 
formations.  The  use  of  only  one-half  the  grid  points,  for  instance, 
during  most  of  the  time  marching  procedure  would  significantly  reduce 
the  computer  time  required  for  an  equivalent  numerical  solution. 


SECTION  VI 


CONCLUSIO  IS  AND  RECOMMENDATIONS 

Numerical  solutions  have  been  obtained  for  two-dimensional  incom¬ 
pressible  turbulent  flow  over  airfoils  near  stall.  The  time  dependent 
Reynolds  averaged  incompressible  Navier-Stokes  equations  in  the  primitive 
variables  of  velocity  and  pressure  together  with  a  Poisson  pressure 
equation  are  numerically  solved.  An  algebraic  eddy  viscosity  approach 
is  modified  for  separated  adverse  pressure  gradient  flows  and  used  to 
model  turbulent  closure  of  the  laminar  separation  bubble  and  the  subsequent 
turbulent  boundary  layer.  A  deficiency  in  the  standard  model  is  detected 
and  corrected  by  using  a  "limiting"  technique.  A  body-fitted  coordinate 
system  is  numerically  transformed  to  a  rectangular  grid  in  the  computational 
plane.  The  set  of  transformed  partial  differential  equations  is  solved 
with  an  implicit  finite  difference  method.  Successive-over-relaxation 
iteration  is  us°d  to  solve  the  system  of  linearized  difference  equations 
at  each  time  step. 

Numerical  solutions  are  presented  for  a  NACA  0012  airfoil  near  stall 
at  a  chord  Reynolds  number  of  170,000.  A  short  laminar  separation  bubble 
located  near  the  upper  surface  suction  pressure  peak  is  obtained.  Computed 
laminar  separation  bubble  characteristics  including  the  criteria  for  sep¬ 
aration,  bubble  type,  and  turbulent  transition  agree  with  empirical  results. 
Surface  mean  pressure  distributions  are  presented  and  found  to  compare 
favorably  with  experimental  data.  The  separation  bubble  is  observed  for 
angles  of  attack  of  5°  and  7.5°.  For  larger  angles  of  attack,  the  small 
bubble  essentially  disappeared  within  the  numerical  resolution  of  the 
streamwise  grid  spacing.  The  steep  leading  edge  suction  pressure  peak  is 
well  defined  for  each  angle  of  attack.  Velocity  profiles  at  four  stations 


4 


S3 


along  the  upper  airfoil  surface  are  compared  with  experimental  results. 

The  experimental  data  were  obtained  at  similar  grid  point  locations 
so  that  interpolation  was  not  required.  The  calculated  lift  and  drag 
coefficients  are  in  excellent  agreement  with  the  experimental  data. 

The  lift  coefficients  are  within  57  of  the  experimental  values  near  stall, 
and  the  computed  drag  coefficients  are  within  10  drag  counts  in  the  region 
of  maximum  lift  to  drag  ratio.  The  effects  of  viscous  separation  on  the 
lift  coefficient  curve  which  produces  a  maximum  value  are  seen  and  con¬ 
trasted  with  a  numerically  obtained  inviscid  potential  solution  with  Kutta 
condition.  The  observed  phenomena  of  trailing  edge  stall  is  predicted 
where  the  rear  separation  point  moves  forward  with  increasing  angle  of 
attack . 

The  sensitivity  of  the  numerical  solution  has  been  examined  in 
several  areas.  A  far  field  potential  flow  boundary  condition  which 
modifies  the  freestream  conditions  for  use  at  a  large  but  finite 
distance  from  the  airfoil  was  considered  for  two  outer  boundary  placements. 
The  study  showed  that  the  use  of  the  infinite  freestream  conditions  at  an 
outer  circular  boundary  of  radius  10  chords  produces  negligible  influence 
on  the  near  flow  field  and  calculated  values  for  the  force  coefficients. 

A  coarse  79x23  grid,  compared  with  the  79x44  grid,  was  used  to  evaluate 
truncation  error.  An  analysis  of  convergence  for  successive-over- 
relaxation  iteration  predicts  an  upper  limit  on  the  time  increment  for 
the  implicit  finite  difference  method.  Numerical  experiments  confirmed 
this  upper  bound.  Several  parameters  within  the  turbulence  model  were 
systematically  varied.  The  only  significant  sensitivity  occurred 
near  the  downstream  side  of  the  laminar  separation  bubble.  Small  changes 
in  turbulent  transition  length  and  location  were  found  to  signif  icant ly 
change  the  flow  field.  This  sensitivity  may  he  related  to  incipient 
turbulent  separation  which  results  in  failure  to  close  the  separation 
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bubble.  This  leading  edge  stall  phenomena  is  observed  at  higher 


Reynolds  numbers  for  the  NACA  0012  airfoil. 

The  agreement  between  the  numerical  solu’ ions  and  the  experimental 
data  and  empirical  results  indicates  that  the  eddy  viscosity  approach 
modified  for  separated  adverse  pressure  gradient  flows  adequately 
models  the  turbulence  on  the  upper  airfoil  surface.  The  near  stall 
airfoil  aerodynamic  characteristics  are  consequently  accurately  pre¬ 
dicted  by  the  numerical  method. 

The  results  of  this  investigation  have  suggested  areas  for  further 
research.  A  more  exhaustive  study  of  the  laminar  separation  bubble 
should  be  accomplished  to  better  understand  the  phenomena  of  leading  edge- 
stall.  The  computation  of  a  geometry  with  detailed  experimental  results 
for  this  purpose  is  desirable.  Further  research  in  adapting  the  multi-  1 

grid  technique  would  significantly  increase  the  computational  efficiency 
and  appears  feasible  as  a  result  of  the  coarse  grid  calculation.  Turbulence 
modelling  is  an  area  of  substantial  interest  and  should  be  pursued  in 
conjunction  with  experimental  investigations.  The  improved  accuracy  of 
current  computational  aerodynamics  also  indicates  that  the  errors  in  ex¬ 
perimental  measurements  should  be  analyzed  in  a  more  rigorous  manner  so 
that  comparisons  can  be  better  evaluated.  The  method  should  also  he 
employed  in  time  dependent  problems  to  further  exploit  the  time  marching 
technique . 
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FIGURE  12  PEAN  FLOW  PRESSURE  CONTOURS  («=  5‘) 


FIGURE  13  JEAN  FLOW  PRESSURE  CONTOURS  (>  =  7.5') 
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FIGURE  20  MEAN  FUDW  VELOCITY  VERSUS  SURFACE  NORMAL  DISTANCE  ICNDlIfNSIONAUZED  BY 
COfPWED  BOUNDARY  LAYER  THICKNESS  AT  CHOFD  LOCATIONS  ON  UPPER  SURFACE 
(a  =  5*) 


FIGURE  21  f€AN  FLOW  VELOCITY  VERSUS  SURFACE  NORMAL  DISTANCE  fCNDIfENSIGNALIZED  BY 

computed  boundary  layer  thicwcss  at  chord  locations  on  lfper  surface 
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FIGURE  22  rt/VJ  FLOW  VELOCITY  VERSUS  SURFACE  ftMAL  DISTANCE  NDNDIIENSI0NAL1ZED  BY 
OOFPUTED  BOUNDARY  LAYER  THICKNESS  AT  CHORD  LOCATIONS  ON  UPPER  SURFACE 
(<*  =  9.5‘) 


FIGURE  25  ITAN  FLOW  VELOCITY  mJS  SURFACE  NORMAL  DISTANCE  NONDirtNSlOHALlZED  BY 
amJTED  BOUNDARY  LAYER  THICKNESS  AT  CHORD  LOCATIONS  ON  UPPER  SURFACE 
<»  =  11.5') 
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FIGURE  2A  WAKE  HAN  FLOW  VELOCITY  PROFILES  AT  CHORE  LOCATIONS  WITH 
ORIGIN  ON  CHORDLttE  (■>  =  5") 
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FIGURE  29  REYNOLDS  STOSS  U’V  CONTOURS  ( 
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FIGURE  30  REYNOLDS  STRESS  U'V'  CONTOURS  (a  =  9.5") 
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FIGURE  31  REYNOLDS  STRESS  iTV'  CONTOURS  (■ 
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FIGURE  36  CONTOUR  INTEGRATION  GEQFETKY  FOR  TWO-DIftNSIONAL  FLOW 
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FIGURE  37  mxiftfl  ERROR  OF  SOR  SOLUTION  FOR  U,  P,  AND  Pb  VERSUS 
fUBER  OF  ITERATIONS  WITH  At  =  0.0006  *«=  8” 
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FIGURE  38  MAXIMUM  ERROR  OF  SOR  SOLUTION  FOR  U,  P,  AND  Pb  VERSUS 
TUBER  OF  ITERATIONS  WITH  At  =  0.001  AND  -  =  8' 
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APPENDIX  A 

GRID  TRANSFORMATION'  RELATIONSHIPS 

Several  of  the  definitions  and  coordinate  transformation  relations 
in  the  transformer  plane  are  summarized  below.  The  notation  used  by 
Thames  (44)  and  Hodge  (53)  is  retained.  The  physical  plane  is  the  (:•:,  y) 
coordinate  system  and  the  computational  plane  is  the  (  :,  ’)  coordinate 
system.  The  following  function  definitions  are  used: 

f(x,  v,  t)  -  a  twice  differentiable  scalar  function. 

_F(x,  y,  t)  =  i  F^(x,  y,  t)  +  j  F?(x,  v,  1 1  a  twice  differentiable 
vector  function  where  i  and  j  are  Cartesian  unit  vectors. 


Transformation  Definitions 

1  =  x~  +  v“  (A  .  1 ) 

r  =  x  +  v  ,vr  (A  .  2 ) 

o  1 

\  -  x“-  +  y  .  (A  .  3) 

D1  =  (.ix  rr  -  2?  x  +  l>tT  _ )  (A.  4) 

T\  =  (Ay  ,  r  -  2;  y  -r  +  -,y_)  (A.  5) 

J  =  x  -V  -xv.  (A .  A ) 

r,  tv  t 

n  =  (y  px  -  x  pJ/J  (A.  7) 

T  =  (xr.D,  -  IVD^/J  (A.M 

D orivatives  of  Scalar  Functions  in  Transformed  Plane 

f  =  •f/'-x  =  (v  f  .  -  v  f  )/.]  < a.ot 

x  '  •  • 

f  =  -f/'.-y  =  (x  J,  -  xrf  .)/.!  (A.  in) 

f  =  f/.-t  =  Of/  t)  ,  r  for  fixed  coordinates  (A.  11) 
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Vector  Operators  in  the  Transformed  Plane 
Grad  lent : 

7f  =  ;(yrfr  -  y  fr  )  i  +  (x  i,  -  y.  f  J  j  /J  (A. 12) 

Divergence : 

7  •  F  =  -,yr  (F1)  -  y  :(F] )  +  x  iF.J  -  x,  (Fj  '7.1  (A. 13) 

Cur  l : 

7  x  F  =  k  ’ v  (F0)  -  v  ,(F0)  -  x/F.)  +  (F,)  7.1  (A.  14) 

-  f  1  -  ,  1  y  1  F 

Lap  lac ian : 

7“f  =  (.  ;f  .  .  -  2rf  j,  +  vf_)/j“  +  (7.J2)7  +  ("  /J~)i  .  (A. 15) 


Unit  Normal  and  Tangent  Vectors  in  Phvsical  Plane 


Normal  to  '  line: 

n(L)  _ 

n  -  •/ 


=  (-V  i  +  x.j)/> 


(A. 16) 


Norm.al  to  f  line: 

n<*  *  7  F/  7  T  =  (vr  i  -  xr  j)/> 


(A. 17) 


Tangent  to  v  line  (increasing  r  direction!: 
t_K  -  n  x  k  =  ( x  jd  +  y  rj )  /  v  •, 


(A. 16) 


Tangent  to  7  line  (.increasing  r,  direction): 

(  t)  (  D)  '  '  _ _ 

ty  =  -ny  x  k  =  (x^i  +  yr  j  )  />  a 


( A  .  1  ? ) 


Integrals  in  the  Transformed  Plane 
Area  Integral: 

^J~Rf(x,  y)  dxdy  =  S?*  f  (x(  e. 


)  ,  y  (  F,  r  ))  ,  J  d  dv 


where  R  is  the  region  R  mapped  into  the  (f,  r.)  plane. 


(A. 20) 
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Vector  Line  Integral: 


_F(x,  y)  ds 


r  ))>  yd  f 
P 


(A. 21) 


where  contour  sr is  any  constant  t]  line  in  the  physical  plane,  s  is  arclength 

alone  s  ,  and  n  is  the  value  of  r  for  the  chosen  contour  s  . 

a  y  d 
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APPENDIX  B 


FINITE  DIFFERENCE  APPROXIMATIONS 


This  appendix  presents  the  finite  difference  approxiinat ions  that  are 
used  in  this  investigation.  The  differences  are  formulated  using  a  func¬ 
tion  f(;,  r.)  in  the  computational  or  transformed  plane.  The  :  and 
spacings  are  assumed  constant  with  value  unity.  The  truncation  error  term 
is  given  in  derivative  form  assuming  L  ?  and  ir  are  equal.  The  approxima¬ 
tions  are  second  order  accurate  unless  otherwise  specified.  Time  dif¬ 
ferences  are  expressed  by  it.  The  superscript  n  denotes  the  n*"'1  time- 
interval  and  is  understood  when  omitted.  The  differences  are  given  for  a 
(  ?,  n)  point  location  denoted  by  subscripts  (i,  j)  and  are  understood  when 
omitted.  Space  derivatives  with  respect  to  only  are  presented  because 
the  corresponding  ■-  derivatives  are  identical  with  subscripts  (i,  j)  reversed. 


First  time  derivative,  backward  difference,  first  order: 
ft  =  (fn  -  fn_1)/it  -  it  ftt/2  +  . .. 

First  derivative,  backward  difference,  first  order: 

fr,  =  <fj  -  fj-l>  "  W2  + 

First  derivative,  central  difference: 


=  (f 


j+1 


-  f 


j-l 


)/2  -  f  / 6  + 


First  derivative,  forward  difference: 


fH  =  (_f j+2  +  4fj+l  -  3V  -  W3  + 


(B.l) 


(E.2) 


( B  .  3) 


(B.i) 


First  derivative,  backward  difference: 


f 


rj 


4f .  .  +  3f ,)  +  f  / 3  +  .  . . 
J-l  J  DV, 


(B .  3) 
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Second  derivative,  central  difference: 


r 


f  =  (f  -  2f .  +  f .  ,)  -  f  / 12  +  .  .  . 

nn  j+i  j  j-i  r:  nn 


Second  derivative,  forward  difference: 

f  =  (-f . . _  +  4f  -  5f.  .  +  2f . )  +  Ilf  / 1 2  + 

CC  J  +  3  1  +  2  1+1  1  rr;t,  . 


Second  derivative,  backward  difference: 


f  =  (f.  _  -  4f.  .  +  5f.  ,  -  2f . )  -  Ilf  / 12  + 
rr,  1-3  i-2  i-l  j'  rrr- 


Cross  derivative,  central  difference: 


f  £h  (fi+l,  j  +  1  f i+1 ,  j-1  +  f i-l ,  j-1  ~  f i-l ,  j+l'/4 


(f?«r,  +  £t™>'24  +  ■■ 


Cross  derivative,  central  in  r  and  forward  in  C  differences: 


f.fr,  =  ;-3(fi+l  ~  £i-l)  +4(fi+l,  j  +  1  -  fi-l,  j+1> 


(fi+l,  j  +  2  f  i-l ,  i+2^  /4  +  (fff 


fr.  +  f  ^r)/6  + 


Cross  derivative,  central  in  C  and  backward  in  n  differences: 
f  '3(fi+l  "  fi-l)  “  4(fi+l,  j-1  "  f  i-l ,  j-P 

+  (fi+i,  j-2  -  fi-i,  j-2)}/4  -  (fm  +  f^)/6  + 


(B.  6) 


(B.7) 


(B .  8) 


(B.9) 

\ 


(B.10) 


(B.ll) 
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APPENDIX  C 


GRID  MODIFICATION  FOR  WAKE  RESOLUTION 

This  appendix  describes  the  technique  that  is  used  to  insert  addi¬ 
tional  grid  lines  in  the  region  of  the  airfoil  wake.  The  body-fitted  coor¬ 
dinate  grid  system  becomes  coarse  at  distances  far  from  the  body.  The 
coarse  grid  may  reduce  the  resolution  of  the  velocity  defect  in  the  viscous 
far  wake  region.  In  order  to  increase  wake  resolution,  eight  additional  ’ 
lines  are  placed  downstream  of  the  airfoil  with  body  points  on  the  rounded 
trailing  edge.  The  original  numerically  generated  grid  has  71  r  and  44  r 
lines.  Linear  interpolation  between  adjoining  ?  lines  of  the  form  x  (I,J) 

=  x ( I — 1 ,  J)  +  x(I,  J)}/2  is  used  to  locate  a  F  line  between  lines  1  =  1 
and  I  =  2,  I  =  2  and  I  =  3,  I  =  69  and  I  =  70,  and  I  =  70  and  I  =  71.  The 
procedure  is  repeated  using  the  new  lines  as  well  except  only  one  additional 
line  is  located  below  the  1=1  cut.  In  this  way  the  fine  grid  spacing  has 
an  asymmetry  for  better  wake  resolution  at  angle  of  attack.  The  final  79 
by  44  grid  is  shown  in  Figure  3.  Table  C.I  gives  the  numbered  1  index 
(F  line)  designation  for  the  grid  systems  in  the  wake  region. 
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APPENDIX  D 


CONTOUR  INTEGRAL  METHOD  FOR  DETERMINING  BODY  FORCE  COEFFICIENTS 

The  determination  of  the  resultant  force  that  a  flow  field  exerts  on 

a  body  usually  involves  the  calculation  of  surface  pressure  and  viscous 

stresses  and  the  summation  of  these  forces  over  the  body  surface.  An 
alternate  approach  is  to  apply  a  control  volume  analysis  to  a  region 
enclosing  the  body.  In  the  foregoing  discussion,  a  control  volume  analysis 
for  a  body  in  a  two  dimensional  flow  is  presented. 

Consider  a  general  body  immersed  in  a  flow  field  w’ith  surface  tg. 
Define  a  fixed  control  volume  V  with  outer  surface  and  inner  surface  0  . 

Conservation  of  linear  momentum  for  the  fluid  in  V  is  expressed  by  the 

Cauchy  Integral  Equation  of  Motion 


(D.l) 


where  .  is  the  density,  v  velocity  field,  b_  body  force  per  unit  mass,  n 

unit  outwara  surface  normal  vector,  T  traction  stress,  and  0  is  the  total 

surface  n  +  o^.  Apply  the  Reynolds  Transport  Theorem  to  the  material 

derivative  term  and  the  Divergence  Theorem  to  the  surface  integral  and  obtain 

t  h 

from  Equation  D.l  for  the  i  Cartesian  vector  component 

I  v—  (,  v .  )  dV  +  I  b  .  (•  v  .  v .  -  T  .  . )  dV  =  0  ( D  •  2 ) 

J  't  l'  J  j  l  J  ji 

V  V 

Apply  the  Divergence  Theorem  to  the  second  volume  integral  and  express  the 
resulting  surface  integral  as  separate  integrals  over  ■“  and  to  obtain 
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/x—  (?  v  . )  dV  +  I  n  .  (p  v  .  v  .  -  T  .  . )  dc  +  j  n  .  (.  v ..  v  .  -  T  .  . )  d  ' 

"t  1  J  J1J  JX  w/  j  '  1  j  JI 


(D.3) 


Now,  on  the  solid  non-porous  bodv  surface  c„ ,  n .  v  =0.  The  integral 

B  j  j 

of  the  resultant  traction  force  over  c  is  the  net  force  exerted  bv  the 

D 

body  on  the  control  volume  fluid.  Therefore,  the  force  _F  exerted  by  the 
fluid  on  the  body  is  given  by 


=  -  S n j 


T.  .  dc 
JX 


Substitute  Equation  D.4  into  D.3  and  use  the  above  surface  condition  to 
obtain  a  general  expression  for  the  force  exerted  by  the  fluid  on  the  body. 


' .  =  -  j  P  v  .  dV  +  /  n  .  (T  .  .  -  o  v .  v  . ' 

l  .-it  J  l  J  J  Ji  xj 


Now  consider  the  case  of  two-dimensional  flow  in  the  x-y  plane.  Then 
Equation  0.5  becomes  the  following  expression  for  the  force  per  length  of 


span  f_ 


'  .  -  -  r—  J  . ' v .  dA  +  /  n  .  (T . .  -  .  v . 

l  t  J  l  J  j  ji  l 


v  . )  ds 
J 


where  A  is  the  cross  sectional  area  of  the  control  volume  and  s,  is  the  outt 

perimeter  of  A  as  shown  in  Figure  36. 

Substitute  for  the  traction  stress  1..  =  -pr  .  .  +  in  Equation  I'.t- 

Ji  J i  J  i 


and  obtain 


.  =  /  n . (-p: . .  +  t . .  -  v.v.lds  -  —  /  .  v .  dA 

1  J  J  Ji  JX  l  J  - 1  J  l 


where  p  is  the  pressure,  •  .  is  the  laminar  viscous  stress  and 

is  the  Kronecker  Delta  Function.  Next  express  the  1  low  variables  in 
terms  of  turbulent  mean  and  fluctuating  variables,  assume  incompress ibl< 
flow,  and  Reynolds  average  Equation  D.7  to  obtain 


1  r  j 


- p  .  .  +  '  .  .  -  ( v . v .  +  v ! v ! )  ds 

Ji  ji  1  J  1  J 


-4,  f  S  dA 


The  term  -  v !  v !  is  the  turbulent  Kevnolds  stress  ...  let  tin  tot.:',  st  re 
t  J  '  t  j  l 

be  defined  as  : .  .  =  '  .  .  +  \  .  .  .  Then  Equation  P .  b  becor.es 
I  J  i  J  i  t  j  i 


r.  /• 

l  J  1 


l  n  .  ■  -p  '  .  .  +  .  .  -  .  v  .  v  .  us 

J  J  J  i  1 .1  l  i  .1 


Introduce  the  following  nond imensional  variable?  in  Equation 


P  -  Pr 


v  . 

\ 

vi =  r~ 


and  t  =  —  ,  and  not*  t  '..it 


1 j i  Tji 


where  r  is  the  freest  rear,  velocitv  and  !  is  a  r*  f  *.  r< 


length  in  the  x  direction.  Also,  define  the  force  coefficient  (  =  'i 

i  I'd 


and  obtain  from  Equation  D.Q 

f ,  =2  S  n  .  -  p 1  .  .  +  , 


,  =2  J  n .  -p '  .  .  +  ,t  .  v .  v .  ds  - 

f  .  s,  J  J  1  '  ,1-  Iji  1  J 


Jv 


where  all  variables  arc  understood  to  be  mean  dimensionless  vari; 

the  turbulent  eddy  viscosity  concept  is  used  where 

"... .  .  =  (■•  +  •  *  C-.v.  +  •  .  v .  ) 

1  J  1  >•  i  J  J  1 


and  define  Ke^  =  +- 


which  transforms  Ecuation  P.  IP  to  the  'cl  levin. 


r,  =2 


-P‘,.  +  o  (  . v .  +  • . v  )  -  v  v .  ds  -  2  - 


,  /v 


<  P .  1!  ) 


Expand  Equation  I).  11  into  x  and  v  component  equations 


I 


/ 


1 

+  Re. 


+ 

1  'X 


!  -  (u  n, 

.•v  .-x  1 


+  uvn„)  ds 


C  . 
i 

v 


(D. 12a) 


(uvr.j  +  vr. .,)  ds 


(  D .  1 2  b  ) 


vliuro  n.  and  n,  ari  the  x  and  y  components  of  the  unit  outward  normal 
vector  n  for  the  outer  boundary  st  of  area  A.  Equations  (I>.  12a)  and  1 

are  expressions  for  tin.'  x  and  v  force  coefficients  exerted  by  a  two- 
dimensional  incompressible  turbulent  flow  on  an  arbitrary  two-dimensional 
body  where  path  s,  encloses  the  body  and  an  eddy  v.scosity  method  is  used 
to  model  turbulence. 


Lift  and  drag  coefficients  are  then  obtained  from  the  force  coefficient 
as  follows  where  ■  is  the  geometric  angle  of  attack. 
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APPENDIX  K 


TIME_1NCREMENT  INFLl'ENCr  ON  SOK  CONVERGENCE 

Implicit  SOK  methods  have  been  shown  by  linear  stability  analyses 
to  be  unconditionally  stable.  However ,  convergence  has  required  a  small 
time  step  size  in  this  investigation  and  elsewhere.  The  following  analy¬ 
sis  indicates  the  relationship  between  the  time  step  size  and  convergence. 

Consider  the  system  of  difference  equations  to  be  written  in  thi 
following  matrix  form: 


AT  =  b  ( L . 1  I 

where  A  is  the  coefficient  matrix,  T  is  the  vector  of  unknown  primitive 
variables  at  time  step  n  with  u,  v,  p  grouped  together  for  each  (i,p 
point  location,  and  b  is  a  vector  of  "constants".  For  example,  the  x  cor 
ponent  of  the  momentum  equation.  Equation  92,  at  the  (i,j)  point  is 

Q  u(5)  =  un .  ^  +  i-t  [  -YETA  (p  .  -  p  .)  +  YXI  (p.  -  p.  .  ) 

lj  ij  i  +  l  j  *  l-l  .)  1  l  j  +  1  1  l  j-1 


+  TC.  (4u 

1C1 

.  -  u 
J 

ic: 

UIV2  + 

w 

(du. 

TV  1 

+  GAMA  (u. 

j+1 

+  u  . 
1 

j-1 

u .  , 
l-l  j+1 

-  u 

(s-1) 

ij 

Pi 

.  )  +  VC  (du  . 


i  TCI 

1 


Ui  tcju  ’T  (4uivi  i 


U i  JV2}  +  RET  ALKA  ("i+l  j  +  "i-1  j ) 


i+l  i+l 


ui+l  i-1  +  “i-1  i-1 


O'.:) 


where  all  coefficients  are  evaluated  at  point  (i,j)  and  primitive  variah 
are  at  time  step  n  unless  given  otherwise,  and  Q  is  givtr.  Ly 
Q  =  1  +  3  .‘.tt  t'C 

Rearrange  Equation  E.2  and  obtain 
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A  sir.ilar  y  component  of  momentum,  equation  is  obtained  if  u  is  replaced  by 
v;  -  YETA  by  XFTA;  and  YXI  bv  -  XXI. 

The  finite  difference  pressure  Equation  94  becomes 
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where  B=2  (ALFA+GAMA)  +  1  (  VV.  +  .VY.)  and  RHS  =  -  ~~  +  [  (EX)  ~+2  (VX)  (FY)  +  (VY ') "  : 

Convergence  theorems  (SI)  for  systems  of  linear  equations  with  constant 
coefficients  show  that  a  necessary  and  sufficient  condition  for  convergence 
requires  that  the  coefficient  matrix  be  diagonally  dominant.  If  this  criteria 
is  applied  to  the  above  quasi-1  inear  system,  the  following  convergence  criteria, 
inequality  E.6  from  the  momentum  equation  and  inequality  E.b  iron  the  pressure 


equation,  are  obtained: 
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The  convergence  criteria,  Equation  E.6,  provides  a  liir.it  on  At  as-  a 
function  of  both  the  grid  geometry  and  flow  field  solution.  The  right  hand 
side  of  the  inequality  has  its  smallest  magnitude  near  the  trailing  edge 
and  increases  rapidly  away  from  the  airfoil.  An  order  of  magnitude  analysis 
of  the  grid  coefficients  and  an  examination  of  the  computational  solutions 
indicate  that  diagonal  dominance  occurs  if  At  <  0.0003  for  points  near  the 
trailing  edge,  t  0.001  elsewhere  near  the  body,  and  At  •  1  in  the  far 

field.  The  pressure  equation  convergence  criteria,  Equation  E.8,  has  no  At 
dependence  and  is  approximately  satisfied  only  at  points  far  from,  the 
body  where  the  left  hand  side  approaches  zero.  These  criteria  can  only 
indicate  a  possible  time  step  restriction  since  the  system  of  equations  has 
variable  coefficients  and  lower  order  nonlinearities. 

A  numerical  experiment  was  conducted  to  determine  the  time  step  effect 

on  convergence.  A  computation  for  the  flow  over  the  NAC.A  0012  airfoil  at 

an  angle  of  attack  of  eight  degrees  was  used.  The  solution  was  advanced 

one  time  step  with  the  time  step  size  and  number  of  iterations  per  time 

step  varied.  The  solution  converged  for  time  steps  of  0.0005  and  0.001  while 

diverging  for  a  tine  step  of  0.005.  The  relative  maximum  error  magnitudes 

for  the  u  component  of  velocity,  surface  pressure  p^ ,  and  field  pressure 

p  as  functions  of  number  of  iterations  are  shown  for  each  tine  step  in 

Figures  37,  38,  and  39.  The  relative  error  is  defined  as  (f  -  f  ,3/f  . 

n  n-1  n 

The  v  component  of  velocity  errors  behave  similarlv  to  the  u  component. 

Next,  the  pressure  was  held  fixed  and  the  convergence  of  the  velocity 
was  monitored  (Figures  37  ,  38,  and  393.  The  rapid  rate  of  convergence  i  s 


observed.  The  u  and  v  velocity  components  were  also  held  fixed  i 
turn  to  test  convergence.  In  both  cases,  convergence  was  similar 
the  general  solution  case,  These  numerical  experiments  indicate  t 
the  convergence  of  pressure  is  slow  and  that  convergence  is  depen 
on  time  step  sizes  similar  to  those  predicted  by  the  simple  an.aly 
above . 


APPENDIX  F 

FAR  FIELD  BOUNDARY  CONDITIONS 

The  far  field  boundary  for  a  numerical  solution  must  usually 
be  at  a  finite  distance  from  the  body  of  interest.  This  constraint 
may  pose  difficulties  if  the  only  known  far  field  boundary  conditions 
are  those  at  infinity.  This  analysis  uses  complex  variable  methods 
for  incompressible  potential  flow  and  develops  approximate  far  field 
boundary  conditions  for  use  at  a  finite  distance. 

Consider  a  circular  cylinder  radius  a  and  center  at  the  origin ,  wit':, 
positive  clockwise  circulation  in  a  uniform  stream  (velocity  Y), 

at  angle  of  attack  _t.  This  coordinate  system  is  identical  with  the 
physical  plane  of  the  airfoil.  The  complex  velocity  is 
..  _ .  (-  ~a " )  i.<  , 

W  (Z)  =  l'rj.  e  -  - ~  —  -rmr —  +  i  —  ~  (F .  1) 

i  z 

where  the  first  term  is  the  uniform  stream,  contribution,  the  second 
term  is  a  doublet  of  strength  2~a~U  ,  and  the  third  term,  is  a  vortex  of 
strength  7. 

Next,  use  the  Kutta-.Joukovski  Theorem  for  lift  and  the  definition 
of  the  lift  coefficient  to  ob  t.  g  in 

r  =  (i/2)  cLrtc  (f . 2) 

where  c  is  the  airfoil  chord. 

Nondimensional  ize  Equation  F.l  velocities  with  the  freestream  velocity 
and  lengths  with  the  airfoil  chord  c  and  substitute  Equation  F.2  into  F.l 


and  obtain 


W  (2)  =  e_1 1  -  (2)' 


+  i 


f  1 


(K.3) 


where  W  is  the  nondimensional  complex  velocity  and  Z  is  the  nondinensional 
complex  variable  Z  =  x  +  iy. 

Airfoil  transformations  such  as  the  Joukowski  and  Treffetz  Trans¬ 
formations  transform  a  cylinder  into  an  airfoil  with  a  chord  of  four 
radii.  Thus,  for  a  given  chord  c,  choose  the  cylinder  radius  to  be  1/4 
c  which  gives  a  scaling  factor  ~  =  1/4  in  the  doublet.  Then  the  complex 
velocity  for  the  flow  around  the  scaled  cylinder  becomes 
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Expression  F.4  can  then  be  used  to  approximate  the  far  field  potential 
flow  over  the  airfoil  of  chord  c  with  circulation  7  in  the  same  physical 
plane . 

Let  the  far  field  boundary  be  a  circle  of  radius  r ,  where 


Z  =  rf  e  \  0  _ 


Substitute  for  Z  in  Equation  F.4  and  obtain 

C, 
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Combine  the  exponential  terms  in  F.5  and  use  the  fact  that  V  =  u-iv, 

where  u  and  v  are  the  x  and  v  velocity  components,  to  obtain 
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In  both  Equations  F.6,  the  first  term  is  the  freestream  velocity  contribution, 
the  second  term  is  the  doublet  contribution,  and  the  third  term  is  the 
vortex  contribution.  As  r^.  approaches  infinity,  the  infinity  freestream 
boundary  conditions  for  u  and  v  are  recovered.  Also,  the  doublet  term 
is  of  higher  order  in  1 / r ^  than  the  vortex  term  for  the  case  of  a  lifting 


airfoil  and  mav  be 


d. 


1  30 


mm mmmmi 


The  far  field  pressure  is  obtained  by  applying  Bernoulli's 


equation  for  irrotational  incompressible  flow  and  using  the  nondimensional 
fo 'm  for  pressure.  The  result  for  the  far  field  boundary  is 

Pf  C)  =  1/2  ;i-(uf:  v f  )  (F.7) 

The  far  field  boundary  conditions  of  liquations  F.6  and  F.7  modify  the 
infinity  freestream  conditions  by  incorporating  effects  on  the  flow  of 
body  thickness  and  circulation  at  a  large  but  finite  distance  from 
the  airfoil. 

Transonic  small  disturbance  theory  for  slender  bodies  and  airfoils 
has  similar  expressions  for  the  far  field  conditions.  Klunker  (101)  has 
obtained  an  asymptotic  solution  applicable  at  large  distances  from  thin 
airfoils.  This  form  has  been  used  as  a  numerical  outer  boundary  condition. 
For  the  limiting  case  of  incompressible  two-dimensional  flow,  the  doublet 
and  vortex  strengths  are  compared  below  with  the  general  forms  found  in 
Equation  F.l. 

The  nondimensional  velocity  potential  doublet  in  the  far  field 
from  Klunker  becomes 


;d  =  2-  FTP  2  f F(s)ds  (F’8) 

where  F(s)  is  the  airfoil  thickness  function.  The  doublet  st-ength  is  thu.- 

2 ^j"F(s)ds.  The  MAC A  four  digit  airfoil  thickness  function  (67)  is 

F (s)  =  (.2969  (-)'*  -  .126  -  -  .3516(-)  +  .2S43  (-)  - 

0.2  c  c  c  c 

.1015  (-)  )  (F.9) 

c 

where  c  is  the  airfoil  chord,  f  is  the  maximum  thickness  as  a  fraction 
of  chord  and  s  is  the  chord  distance,  0  <_  s  c . 

Substitute  the  thickness  function  Equation  F.9  into  the  doublet  strength 
expression  and  obtain  the  following  doublet  strength  from  transonic  small 


disturbance  theorv: 


0.685  fc' 


(F. 10) 
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For  the  case  of  a  NACA  0012  airfoil,  hp  =  0.0822c".  The  corresponding 

doublet  strength  in  the  simple  potential  flow  model  from  Equation  F.l  is  2"' 

c  2  2 

(■£-)  or  .3972c  .  Klunker  points  out  that  the  doublet  may  be  neglected  in 

the  far  field  for  finite  circulation  flows  because  the  vortex  contribution 

1  1  2 

is  of  order  (—  )  compared  with  an  order  (—  )  for  the  doublet. 

rf  rf 

The  vortex  potential  term  of  Klunker  for  the  case  of  incompressible 
two-dimensional  flow  is 

$v  =  sgn  (y)  +  tan  1  (~)  ;  (F.ll) 

where  F  is  the  clockwise  circulation,  the  inverse  tangent  function  is 
defined  on  the  interval  (-  —  — )  ,  and  the  angle  within  braces  is  defined 

^  J  ^ 

clockwise  from  the  negative  x  axis.  This  potential  function  is  identical 
with  the  vortex  function  in  Equation  F.l  where  the  clockwise  strength  s  again 
F  and  the  angular  orientation  is  the  right  hand  polar  coordinate  de,  ed 
counterclockwise  from  the  positive  x  axis. 

Thus  the  far  field  potential  for  each  method  differs  only  in  the 
thickness  or  doublet  term.  Furthermore,  for  flows  with  finite  cir--  ration 
this  difference  is  negligible  because  the  vortex  term  dominates. 
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APPENDIX  C 


NUMERICAL  INVISCID  FLOW  RESULTS 


A  numerical  solution  for  the  inviscid  flow  over  the  NACA  0012 
airfoil  is  obtained  for  comparison  with  both  the  Kaiver-Stokes  computation 
and  the  experimental  data.  The  numerical  algorithm,  developed  by  Thames 
(44),  which  solves  Laplace's  equation  for  the  stream  function  in  two- 
dimensional  flow  is  used.  The  method  uses  SOR  iteration  to  solve  the 
system  of  equations  written  in  terms  of  body-fitted  coordinates.  The 
stream  function  value  on  the  body  can  either  be  specified  or  computed  by 
imposing  a  given  circulation  or  a  Kutta  condition. 

The  stream  function  equation  for  two-dimensional  flow  is 


V  +  V'  =0  (G.l) 

xx  yy 

Rewrite  Equation  G.l  in  terms  of  the  transformed  variables  in  the 
computational  plane  using  expressions  in  Appendix  A  to  obtain 


a  v  r  , 


+a  if 


+  OVr  +T  v  ^  =  0 
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The  boundary  condition  for  the  stream  function  on  the  outer  boundary,  whir 
is  a  radius  of  ten  chords,  is  the  freestream  value 


y^  cos  -  Xj  sin  a 


(i  .  :,i> 


where  (x^.,  y^)  are  coordinates  on  the  outer  computational  boundary  a:  ' 
is  the  angle  of  attack.  Hodge  and  Ghia  (82)  have  showmi  by  an 
inviscid  analysis  that  this  boundary  condition  induces  an  error  in  lift 
of  less  than  IX  for  angles  of  attack  less  than  10°.  The  boundary 
condition  on  the  body  is 

C.  =  constant  (G.Tb' 

b 

The  constant  can  be  specified  or  determined  by  imposing  a  given  circula¬ 
tion  or  a  Kutta  condition.  The  Kutta  condition  which  matches  the  upper 
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and  lower  trailing  edge  surface  velocities  by  extrapolation  is  select¬ 


ed.  The  details  for  each  option  are  given  by  Thanes  (44). 

The  derivatives  in  the  transfcrmed  Equation  G.2  are  approximated 
by  central  difference  formulas  found  in  Appendix  E.  The  resulting 
system  of  equations  is  solved  by  SOR  as  described  by  Thames  (44).  S  lo¬ 
tions  are  obtained  for  0°,  5%  7.5°,  9.5%  and  11.5°  for  the  NAC.A  r':  1 J 
airfoil  section.  The  body  surface  pressure  coefficients  can  be  calculated 
from  the  stream  function  solution  as  follows.  The  body  pressure  coefl'icie 
is  defined  as 

C  =  1  -  (u“  +  v'  )  (('.  4) 

P 

where  at  a  location  on  the  body  u  and  v  are  the  x  and  y  components  of  tin 
flow  velocity  nondimensional  ized  by  the  freestrear.  velocity.  l'sL-  tin. 
definition  of  the  stream  function  to  obtain 


Then  rev/rite  the  Equation  G.5  in  terms  of  the  transformed  variables  and 
note  that  1  r  =0  on  the  body  to  give 

C  =  1  -  t.  (’I)''  (G.6) 

p  J  t. 

where  all  quantities  are  evaluated  on  the  body  surface.  The  body  pressure 
coefficients,  calculated  using  Equation  G.6  with  second  order  one-sided 
differences,  are  plotted  in  Figures  16,  17,  18,  19,  and  33.  A  typical 
streamline  plot  for  a  =  11.5°  is  shown  in  Figure  40.  The  numerical  in- 
viscid  solution  is  also  compared  with  an  analytical  solution  (67)  which 
uses  Theodorsen’s  (12)  method.  The  excellent  agreement  is  seen  in 
Figure  40  where  both  solutions  are  plotted. 
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APPENDIX  H 


EXPERIMENTAL  DATA 


This  appendix  presents  the  unpublished  data  which  have  been  obtained 
experimentally  at  the  Air  Force  Frank  J.  Seiler  Research  Laboratory,  ISA:' 
Academy,  Colorado,  by  laboratory  personnel  and  used  in  this  work. 

Airfoil  upper  and  lower  surface  mean  pressure  coefficient  measure¬ 
ments  for  nominal  attack  angles  of  0°,  5°,  8°,  10°,  and  12°  at  a  chord 
Reynolds  number  of  170,000  are  given  in  Table  H.I.  The  pressure  coefficient 
is  defined  as  the  static  pressure,  relative  to  freest ream  static  pressure, 
nondimensional ized  bv  the  freestream  dynamic  pressure  as  follows:  C  = 

'  -  1  rv 


'p  -  p)/.  ■  l'~  .  The  two  sets  of  hot  wire  anenometrv  velocity  field  data 
near  the  upper  and  lower  airfoil  surfaces  on  designated  paths  are  presented 
in  Table  H.1I.  The  velocity  field  and  Reynolds  stress  u'v'  data  for  the 
near  wake  region  are  presented  in  Table  H.1II.  The  data  have  been  non- 
dimensional  i zed  by  the  measured  freestream  velocity.  The  spatial  locations 
in  the  physical  coordinate  system  for  each  path  are  given  in  Table  H . IV 
where  location  (0,0)  is  the  airfoil  center  and  the  x  axis  lies  along  the 
airfoil  chord. 

The  spatial  measurement  locations  correspond  with  computational  grid 
point  locations  on  constant  r  lines  or  '  lines  in  the  physical  plane.  The 
experimental  lead  time  necessitated  the  use  of  the  original  71  x  44  grid 
system.  The  I,  J  index  notation  designates,  however,  the  point  indexing 
(  r,r)  in  the  final  79  x  44  grid  described  in  Appendix  C. 
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TABLE  H.I  Experimental  Surface  Pressure  Coefficients  for  NACA  0012 
At  Various  Angles  of  Attack,  Re  =  170,000 
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pressure  are  used.  Turbulence  is  modelled  with  an  algebraic  eddy 
viscosity  technique  modified  for  separated  adverse  pressure  gradient 
flows.  The  set  of  transformed  partial  differential  equations  is.  solved 
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with  an  implicit  finite  difference  method.  Numerical  solutions  for 
NACA  0012  airfoil  near  stall  at  a  chord  Reynolds  number  of  170,000 
are  favorably  compared  with  surface  pressure  and  velocity  field 
measurements.  A  small  laminar  separation  bubble  near  the  suction 
pressure  peak  is  observed.  Computed  lift  and  drag  coefficients 
agree  well  with  experimental  values. 
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